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Abstract 

The  interaction  of  submerged  manuevering  vehicles  with  the  free  surface  is  a  topic  of  im- 
portance in  ocean  engineering.  Present  methods  for  estimating  the  forces  and  motions  due 
to  wave-body  interactions  are  limited  in  their  ability  to  efficiently  compute  forces  and  mo- 
tions of  interest  when  free  surface  slopes  are  steep  or  body  motions  are  large.  Methods  are 
required  which  contribute  to  overcoming  these  difficulties. 

This  thesis  presents  a  method  for  analyzing  the  nonlinear  interaction  of  non-breaking 
waves  with  submerged  bodies  of  arbitrary  geometries  undergoing  arbitrary  motions.  The 
method  couples  a  high  order  spectral  representation  of  the  free  surface  profile  and  free 
surface  potential  with  a  boundary  element  representation  of  the  body  to  effect  a  time  domain 
solution  for  specified  initial-boundary  value  problems.  The  spectral  representation  of  free 
surface  quantities  facilitates  the  use  of  fast  transform  techniques  for  rapidly  computing 
free  surface  quantities.  The  boundary  element  repesentation  of  the  body  enables  bodies  of 
arbitrary  geometry  to  be  modelled.  Fully  nonlinear  free  surface  boundary  conditions  are 
used  as  time  evolution  equations  for  the  free  surface  profile  and  potential. 

The  effectiveness  of  the  method  is  demonstrated  through  the  solution  of  three  classes 
of  two-dimensional  problems:  (1)  diffraction  of  incident  waves  by  a  stationary  cylinder,  (2) 
wave  radiation  by  a  cylinder  undergoing  forced  oscillations,  and  (3)  combined  diffraction 
and  radiation  by  a  submerged  neutrally  buoyant  cylinder  free  to  respond  to  incident  waves. 
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Chapter  1 


Introduction  and  Background 


This  thesis  presents  a  method  for  solving  wave  body  interaction  problems.  The  method 
developed  is  appropriate  to  the  analysis  of  the  nonlinear  interaction  of  non-breaking  waves 
with  fully  submerged  bodies  of  arbitrary  geometry  undergoing  arbitrary  motions.  The 
method  is  suitable  for  the  study  of  three  or  two  dimensional  problems.  The  efficacy  of  the 
method  is  demonstrated  through  the  study  of  three  classes  of  two  dimensional  problems: 

•  The  diffraction  of  an  incident  Stokes'  wave  by  a  fully  submerged  cylinder  of  arbitrary 
shape; 

•  The  wave  radiation  by  a  fully  submerged  circular  cylinder  undergoing  large  amplitude 
forced  oscillatory  motion;  and, 

•  The  response  of  a  neutrally  buoyant  circular  cylinder  to  an  incident  Stokes'  wave 
(combined  radiation  and  diffraction). 

The  purpose  of  this  chapter  is  to  explain  the  motivation,  background,  and  format  for 
this  thesis. 
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1.1  Motivation 

This  thesis  is  motivated  by  the  desire  to  efficiently  estimate  and  understand  the  interaction 
of  fully  submerged  bodies  with  the  free  surface.  A  time  domain  solution  for  the  forces 
resulting  from  the  interaction  of  the  body  with  waves  is  necessary  to  perform  maneuver- 
ing simulations  of  bodies  operating  near  the  free  surface.  Present  methods  for  estimating 
forces  due  to  wave-body  interactions  limit  the  computing  speed  of  maneuvering  simulation 
calculations  or  include  linearizing  assumptions  which  limit  the  ability  to  predict  forces  and 
motions  of  interest  when  free  surface  slopes  are  steep  or  unsteady  body  motions  are  large. 
This  thesis  presents  a  method  which  overcomes  some  of  these  shortcomings. 

1.2  Background 

The  Physical  Problem 

The  complete  physical  problem  for  a  body  operating  in  proximity  to  a  free  surface  is 
dominated  by  inertial,  viscous,  and  gravitational  forces  [33].  Other  mechanisms,  such  as 
surface  tension,  exist  but  are  insignificant.  A  mathematical  representation  of  the  complete 
physical  problem  amenable  to  efficient  solution  does  not  exist.  The  problem  includes  the 
possibility  of  separated  flows  and  breaking  waves,  phenomena  difficult  to  analyze.  Numerical 
methods  enable  a  range  of  problems  that  include  some  of  these  phenomema  to  be  solved, 
however  the  computing  time  required  by  some  methods  to  solve  these  problems  for  arbitrary 
bodies  is  prohibitively  excessive.  For  example,  the  solution  of  the  Navier-Stokes  equations 
for  the  flow  around  a  submarine  in  an  unbounded  fluid  for  six  body  lengths  of  travel  requires 
as  much  as  100  hours  of  CPU  time  on  an  IBM  Model  590  workstation  with  512  MB  of  RAM 
[27].    Numerical  methods  are  needed  which  can  characterize  the  relevant  fluid  and  body 
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behavior  with  improved  computational  efficiency. 

The  relative  importance  of  the  principal  force  mechanisms  varies  as  the  parameters 
characterizing  the  problem  change.  Dimensional  analysis  is  used  to  define  the  relationships 
between  the  force  mechanisms  and  identify  ranges  of  problem  parameters  over  which  one  or 
more  of  the  force  mechanisms  may  be  ignored.  For  a  typical  marine  vessel  moving  beneath 
the  free  surface  the  ratio  of  inertial  to  viscous  forces  is  large  outside  the  boundary  layer. 
The  fluid  can  be  treated  as  inviscid  and  the  governing  mathematical  representation  of  the 
problem  appropriately  simplified  [33].  If  the  problem  is  further  restricted  to  considering 
only  non-breaking  waves,  a  wide  range  of  practical  problems  can  be  solved.  In  the  following 
section,  the  existing  methods  for  solving  wave-body  interaction  problems  for  inviscid  fluids 
under  non-breaking  waves  are  briefly  discussed. 

Solution  Methods 

The  simplest  solution  method  is  that  proposed  by  Faltinsen  [9].  This  method  approxi- 
mates the  hydrodynamic  force  as  the  sum  of  (1)  the  Froude-Krylov  Force  (Ffk),  the  force 
due  to  the  unsteady  pressure  field  associated  with  the  undisturbed  incident  wave  potential; 
(2)  a  force  proportional  to  the  product  of  the  body's  added  mass  and  the  acceleration  vector 
of  the  undisturbed  wave  field;  and  (3)  the  added  mass  force,  the  force  resulting  from  the 
body's  acceleration.  This  approach  has  been  used  [5]  for  simulating  the  forces  associated 
with  undersea  vehicles  operating  near  the  free  surface.  The  shortcomings  of  this  type  of 
approach  are  discussed  by  Newman  [32]  and  Milgram  [29].  This  approach  can  only  be 
accurate  if 

•  The  Froude  Numbers  based  on  submergence  depth,  Fh,  and  body  length,  Fi  are 
small: 
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Fh  =  ^=    «1    and    Fl  =  ~j=    «1  (1.1) 


where  Vs  is  the  body  velocity,  g  is  the  gravitational  constant,  H  is  the  body's  mean 
submergence,  and  L  is  the  body's  length;  and, 

•  The  body  is  small  compared  to  the  wavelength. 

Because  of  these  restrictions,  the  method  is  not  suitable  for  many  combinations  of  marine 
vessels  and  operating  scenarios  of  interest. 

Numerical  solution  methods  for  wave-body  interaction  problems  exist  in  both  the  fre- 
quency and  time  domains.  The  application  of  the  principle  of  linear  superposition  to  ship 
motions  in  irregular  seas  was  motivated  by  the  work  of  John  [15]  and  St.  Denis  and  Pierson 
[44].  The  former,  as  noted  by  Korsmeyer  [19],  enabled  the  problem  to  be  decomposed  into 
separate  radiation  and  diffraction  problems,  and  the  latter,  as  noted  by  Salvesen,  et  al. 
[40],  enabled  ship  responses  in  irrregular  seas  to  be  evaluated  as  the  sum  of  the  responses  at 
all  frequencies  of  regular  waves.  Jointly  these  works  motivated  the  development  of  a  wide 
body  of  linearized  ship  motion  analysis  methods  in  the  frequency  domain  using  strip  theory 
and  slender  body  theory  to  characterize  the  wave-body  interaction.  The  work  of  Salvesen, 
Tuck,  and  Faltinsen  [40]  is  a  thorough  application  of  strip  theory  to  ship  motion  analysis. 
Ogilvie  [36]  provides  a  thorough  development  of  slender  body  theory  and  its  application  to 
ship  motions. 

Improvements  in  digital  processing  capabilities  enabled  the  development  of  three  dimen- 
sional analysis  methods  for  arbitrary  wave  induced  ship  motions  [2].  Solution  methods  have 
evolved  systematically,  from  the  solution  to  the  linearized  problem  for  stationary  bodies  or 
bodies  undergoing  forced  motion  to  the  solution  for  bodies  undergoing  nonlinear  unsteady 

16 


motions  with  forward  speed.  Both  frequency  and  time  domain  solution  methods  have  de- 
veloped for  such  problems.  The  solution  to  the  specified  boundary  value  or  initial  value 
problems  is  found  by  either  volume  methods  (finite  difference  and  finite  element  methods) 
or  boundary  integral  equation  methods  (boundary  element  or  "panel"  methods)  [7].  The 
present  preferred  methods  in  ship  hydrodynamics  are  the  boundary  integral  equation  meth- 
ods (BIEM).  Kring  [20]  summarizes  the  development  of  three  dimensional  panel  methods, 
with  emphasis  in  his  summary  on  the  distinction  between  methods  which  use  the  transient 
free  surface  Green  function  (  Korsmeyer  [19]  and  Bingham  [2])  and  those  which  use  the 
Rankine  source  (Nakos  [31]  and  Kring  [20]).  Time  domain  simulations  of  wave-induced  mo- 
tions are  of  immediate  interest  for  their  applicability  to  vehicle  motion  control  and  control 
system  design. 

Dommermuth  [7]  and  Liu  [22]  developed  a  high-order  spectral  method  (HOSM)  for  the 
efficient  time-domain  analysis  of  non-linear  wave-body  interactions.  The  high-order  spectral 
method,  through  the  use  of  Fast  Fourier  Transform  (FFT)  techniques,  greatly  improves  the 
computational  efficiency  of  the  numerical  analysis  of  wave-body  interactions. 

The  solution  method  developed  in  this  thesis  couples  the  boundary  integral  representa- 
tion of  the  body  with  a  spectral  representation  of  free  surface  quantities.  To  demonstrate 
the  efficacy  of  the  method,  the  three  classes  of  problem  identifed  previously,  radiation, 
diffraction,  and  combined  radiation  and  diffraction,  are  studied  for  two  dimensional  bodies. 
The  solutions  are  compared  with  solutions  generated  by  different  methods  for  similar  two 
dimensional  problems.  These  other  methods  will  be  briefly  summarized  here.  A  similar 
summary  of  some  of  these  methods  and  related  methods  can  be  found  in  Dommermuth  [7] 
and  Liu  [22].  More  detailed  discussions  of  the  methods  and  the  solutions  generated  will  be 
discussed  in  subsequent  chapters  as  appropriate. 
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Frequency  domain  solutions  to  wave  diffraction  and  radiation  by  circular  cylinders  were 
developed  by  Ursell  [47]  and  Ogilvie  [35].  Ursell  used  a  multipole  expansion  of  the  pulsating 
source  potential  and  a  series  expansion,  in  powers  of  the  wave  slope,  for  the  coefficients 
in  the  multipole  expansion,  to  derive  a  system  of  equations  to  uniquely  determine  the 
linear  solution  to  the  interaction  between  waves  and  a  fully  submerged  stationary  circular 
cylinder.  Ogilvie  extended  Ursell's  work  to  compute  the  first-order  unsteady  force  and  the 
second-order  steady  force  for  the  three  problems  of  interest  in  this  thesis.  Wu  [53]  extended 
Ogilvie's  work  using  the  multipole  expansion  and  the  free  surface  boundary  condition  for  the 
second  order  potential  to  demonstrate  that  to  second  order,  a  circular  cylinder  undergoing 
a  circular  orbit  generates  waves  in  one  direction  only.  Ogilvie  previously  demonstrated  this 
result  at  first  order.  Using  a  linear  free  surface  boundary  condition  and  the  exact  body 
boundary  condition,  Wu  [52]  found  that  a  circular  cylinder  undergoing  large-amplitude 
circular  motion  can  generate  waves  travelling  in  both  directions  though  the  waves  travelling 
in  the  "upstream"  direction  were  at  most  a  third  order  quantity. 

Chapman  [4]  analyzed  large  amplitude  transient  motions  of  two-dimensional  bodies  with 
a  linearized  free  surface  boundary  condition  and  exact  body  boundary  conditions.  Chapman 
represented  the  wave  field  as  a  finite  sum  of  eigenfunctions  and  the  body  as  a  distribution  of 
sources  and  negative  image  sources  reflected  about  the  mean  free  surface.  For  Chapman's 
method  the  periodic  boundary  conditions  required  by  the  spectral  representation  of  free 
surface  quantities  are  ignored  in  the  representation  of  the  body  potential,  the  harmonic 
components  used  to  represent  the  free  surface  quantities  are  nonuniform  and  can  not  be 
determined  by  FFT,  and  the  method  is  applicable  to  flows  described  by  linear  free  surface 
boundary  conditions. 
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Grue  and  Palm  [11]  studied  the  reflection  of  waves  by  submerged  cylinders  of  circular  and 
elliptic  cross  section  using  boundary  integral  equations.  In  their  formulation  the  solution 
is  expressed  as  a  distribution  of  vortices  on  the  body.  Grue  and  Palm  used  the  linear  free 
surface  boundary  conditions  in  their  formulation  and  confirmed  the  results  of  Ursell  and 
Ogilvie  (no  reflection  from  a  circular  cylinder)  and  computed  the  reflection  power  for  the 
elliptic  contour.  Their  analyses  of  cylinders  of  elliptic  cross  section  are  restricted  to  deeply 
submerged  cylinders. 

Stansby  and  Slaouti  [45]  developed  a  time  domain  solution  method  for  evaluating  non- 
linear wave-body  interactions  with  the  free  surface  represented  as  a  periodic  vortex  sheet 
and  the  body  as  a  distribution  of  periodic  sources.  The  focus  of  their  work  was  on  the  first 
order  oscillatory  forces  and  wave  motions  only  and  they  demonstrated  reasonable  agreement 
with  the  first  order  solution  of  Ogilvie,  in  both  the  amplitude  and  phase  of  the  oscillatory 
forces. 

A  numerical  solution  to  the  second-order  wave  diffraction  problem  in  the  frequency 
domain  was  provided  by  Vada  [48].  Vada  derived  the  boundary  conditions  for  the  first  and 
second  order  diffraction  potentials  using  perturbation  expansions  of  the  boundary  conditions 
and  solved  the  first  and  second  order  problems  using  a  boundary  element  method.  Vada  used 
the  frequency  domain  free  surface  Green  function  provided  by  Wehausen  and  Laitone  [50] 
for  the  first  order  problem  and  a  suitable  modification  of  it  for  the  second  order  problem. 
The  results  generated  by  this  method  compare  favorably  with  the  results  of  Ogilvie  and 
experimental  results  [3].  Vada  did  not  solve  for  the  time-independent  part  of  the  second 
order  potential  and  restricted  his  study  to  stationary  bodies. 

The  Mixed  Eulerian-Lagrangian  (MEL)  method  has  been  applied  to  the  study  of  non- 
linear wave  body  interactions  by  a  number  of  authors.  Vinje  and  Brevig  [49]  and  Cointe  [6] 
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studied  the  interaction  of  cylinders  with  non-linear  free  surface  boundary  conditions.  Vinje 
and  Brevig  demonstrated  the  suitability  of  the  approach  to  the  study  of  wave-body  inter- 
actions. Cointe  studied  the  diffraction  and  radiation  problems  in  detail  and  demonstrated 
the  ability  of  non-linear  solution  methods  to  capture  behavior  not  predicted  by  linear  and 
second  order  theories.  However,  the  Lagrangian  representation  of  free  surface  quantities  is 
not  amenable  to  solution  by  fast  transform  techniques. 

Vinje  and  Brevig's  work,  which  used  Cauchy's  integral  theorem  to  represent  the  velocity 
potential  and  stream  function  ,  was  used  by  Jagannathan  [14]  to  study  radiation  conditions 
for  the  MEL  method.  In  the  process,  Jagannathan  produced  time  simulations  for  the  forces 
on  a  submerged  circular  cylinder.  The  force  results  vary  significantly  from  theortical  and 
experimental  results  for  wave  diffraction  but  show  better  agreement  for  forced  motions.  Ad- 
ditionally, the  difference  between  Jagannathan's  results  and  experimental  results  increases 
as  submergence  decreases. 

Silva  and  Peregrine  [43]  also  used  the  MEL  approach  to  study  nonlinear  wave-body 
interaction  problems  for  submerged  fixed  bodies  and  submerged  bodies  undergoing  forced 
oscillatory  motions.  Silva  and  Peregrine  demonstrated  the  influence  of  finite  depth  on  the 
forces  and  free  surface  profiles  associated  with  forced  body  motion.  Silva  and  Peregrine's 
approach  is  not  extendable  to  three  dimensions. 

The  high-order  spectral  methods  (HOSM)  of  Dommermuth  and  Liu  have  been  used  to 
study  nonlinear  wave-body  interactions  in  both  three  and  two  dimensions.  By  performing 
computations  at  orders  higher  than  second  order,  Dommermuth  and  Liu  were  able  to  com- 
pute a  steady  horizontal  drift  force  acting  on  a  fixed  submerged  circular  cylinder.  This  force 
is  not  predicted  by  second  order  theories  but  has  been  measured  experimentally  [30]  and 
computed  by  fully  nonlinear  methods  [6].  The  HOSM  uses  the  FFT  to  rapidly  compute  the 
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body  and  free  surface  potential.  The  HOSM  requires  that  the  body  be  capable  of  being 
efficiently  represented  by  series  of  orthogonal  functions.  Additionally,  the  HOSM  represents 
the  body  potential  as  a  Taylor  series  expansion  about  the  mean  body  position.  Therefore, 
the  method  can  not  be  used  to  model  a  body  undergoing  large  amplitude  unsteady  motions. 

1.3     Thesis  Format 

This  thesis  consists  of  seven  chapters.  The  mathematical  formulation  of  the  combined 
high-order  spectral  boundary  integral  equation  method  is  provided  as  Chapter  2.  The 
numerical  implementation  of  the  method  is  provided  in  Chapter  3.  The  results  for  the  radi- 
ation, diffraction,  and  combined  problems  are  provided  in  Chapters  4,  5,  and  6  respectively. 
Chapter  7  provides  conclusions  and  recommendations  for  future  work. 
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Chapter  2 


Approach 


This  chapter  presents  the  method  developed  for  evaluating  the  interaction  of  a  fully  sub- 
merged body  with  the  free  surface.  The  problem  statement,  the  solution  formulation,  and 
issues  associated  with  the  formulation  are  presented. 

2.1     Introduction 

The  motivation  for  the  mathematical  formulation  to  follow  is  the  solution  of  the  equation 
of  motion  for  a  submerged  body  subject  to  both  deterministic  and  stochastic  forcing.  The 
deterministic  forcing  is  the  result  of  control  surface  motions  and  operation  of  propulsive 
devices.  The  stochastic  forcing  is  due  to  the  interaction  of  the  body  and  its  appendages 
with  waves.  The  problem  formulation  presented  here  provides  an  efficient  method  for  com- 
puting the  forces  and  moments  on  a  fully  submerged  body  of  arbitrary  geometry  due  to  its 
interaction  with  waves. 

The  formulation  provides  a  mathematical  model  for  the  physical  problem  represented 
by  Figure  2-1. 
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Figure  2-1:  Coordinate  system  definition 


2.2     The  Initial-Boundary  Value  Problem 

Field  Equation 

For  a  large  portion  of  the  range  of  Froude  numbers  and  Reynold's  numbers  at  which  marine 
vessels  operate  while  in  proximity  to  the  free  surface,  viscous  effects  are  confined  to  a 
thin  layer  of  fluid  near  the  vessels'  surfaces.  Prandtl's  boundary  layer  theory  motivates 
and  provides  justification  for  neglecting  viscous  effects  in  the  study  of  selected  problems 
in  hydrodynamics  [33].  Chaplin  [3]  demonstrated  experimentally  that  in  the  absence  of 
viscous  fluid  effects  significant  non-linear  wave  forces  exist. 

Wave-body  interactions  are  dominated  by  fluid  flows  which  can  idealized  as  inviscid  and 
irrotational.  Therefore  a  portion  of  the  non-linear  aspects  of  wave-body  interactions  can  be 
studied  using  potential  flow  theory. 
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For  an  incompressible  fluid  of  constant  density  the  equation  for  conservation  of  mass  is 
the  continuity  equation 

V-V  =  0  (2.1) 

For  an  ideal  fluid  initially  at  rest  and  subject  to  conservative  forces  the  fluid  rotation  rate 
is  zero  and 

V  x  V  =  0  (2.2) 

The  velocity  vector  for  an  ideal  irrotational  flow  can  be  defined  as  the  gradient  of  a 
scalar  potential 

V  =  V$  (2.3) 

Equations  (2.1)  and  (2.3)  are  combined  to  yield  Laplace's  equation,  the  governing  partial 
differential  equation  for  the  velocity  potential: 

V2$  =  0  (2.4) 
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Boundary  Conditions 

The  following  boundary  conditions  are  imposed  on  the  physical  problem  represented  by 
Figure  2-1: 
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on  y  =  n  (2.9) 

-(p-Pa)  =  ^  +  ;:V</)-V0  +  5y  =  O  on  y  =  rj  (2.10) 


where  t/(x,  ^,  t)  is  the  shape  of  the  free  surface,  U  is  the  body  velocity,  and  2L  is  the  domain 
length. 

In  the  absence  of  a  bottom: 

V$  -+  0  as  y  ->  -co  (2.11) 

The  specified  form  of  the  kinematic  and  dynamic  free  surface  boundary  conditions  (equa- 
tions (2.9)  and  (2.10)  respectively)  determines  the  type  of  free  surface  flows  represented. 
An  Eulerian  form  of  the  free  surface  boundary  conditions  capable  of  describing  non-linear 
single-valued  free  surface  flows  is  used  in  this  thesis.  The  chosen  form  was  developed  by 
Zakharov  [55]  for  the  purpose  of  evaluating  the  stability  of  dispersive  surface  waves.   The 
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"Zakharov  Equations"  express  the  free  surface  boundary  conditions  as  functions  of  the 
canonical  variables  $s(x,z,t)  =  $(x,r],  z,t)  and  77,  and  the  vertical  velocity  on  the  free 
surface,  $y(x,r],z,t): 

Vt  =  (l+r]2x)$y(x,T1,t)+$sxr)x  (2.12) 


*f  =  -m  ~  l(*x)2  +  5(1  +  i)«J(«,  V,  t)  (2.13) 


Spatial  periodicity  (equations  (2.7)  and  (2.8))  is  specified.     The  significance  of  this 
condition  is  explained  in  chapter  3. 

Initial  Conditions 

The  initial  conditions  for  the  the  shape  of  the  free  surface,  r}(x,z,0),  and  the  free  surface 
potential,  $s(x,z,Q),  are  specified  as  those  associated  with  Stokes'  waves. 

2.3     Solution  Method 

The  Zakharov  forms  of  the  kinematic  and  dynamic  free  surface  boundary  conditions  are 
integrated  in  time  to  update  the  free  surface  shape  and  free  surface  potential  using  the 
specified  initial  conditions  and  the  solution  of  the  boundary  value  problem.  A  combined 
high-order  spectral  method  (HOSM)  and  boundary  integral  equation  method  (BIEM)  is 
used  to  solve  the  specified  initial-boundary-value  problem.  The  Zakharov  equations  were 
used  in  the  HOSM  developed  by  Dommermuth  and  Yue  [8]  and  West,  et.  al.  [51]  to 
solve  non-linear  wave-wave  interaction  problems  and  by  Liu,  Dommermuth,  and  Yue  [23]  to 
evaluate  non-linear  wave-body  interactions.   Olmez  [37]  used  the  Zakharov  equations  and 
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a  boundary  integral  equation  method  to  analyze  non-linear  wave- wave  interaction  problems. 

The  high  order  spectral  method  (HOSM)  and  (BIEM)  are  summarized  in  the  following 
sections. 


High-Order  Spectral  Method 

The  HOSM  summarized  here  is  that  of  Liu,  Dommermuth,  and  Yue  [23]  for  solving  nonlinear 
wave-body  interaction  problems.  In  the  HOSM,  and  in  the  combined  HOSM-BIEM  solution 
method  developed  herein,  the  velocity  potential  is  expanded  in  a  perturbation  series  up  to 
a  specified  order  M: 

M 

•<m>(*,y,*,t)  =  Y,  *(m)(*.y.*.«)  (2-14) 

m=l 

(  )(m)  denotes  a  quantity  of  order  (e)(m)  where  e  is  a  measure  of  wave  steepness,  kA. 
Each  perturbation  potential,  &m}(x,  y,z,£),  is  expanded  in  a  Taylor  series  about  y  =  0 
such  that: 


M    M-m     I  n   I 
*S{x,Z,t)   =  •(*,!!,*,«)   =:    £      E     ^Sr*(m,(*»tt*.*)  (2-15) 


m=\    1=0 


For  the  specified  initial  conditions,  equation  (2.15)  establishes  a  sequence  of  Dirichlet 
boundary  conditions  for  each  perturbation  potential,  &m)(x,0, z,t): 

${l){x,0,z,t)=<l>s(x,z,t)  (2.16) 
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m—  1    I    <v 

$(m)(x,0,M)  =  -£^^(m~0(*AM),  m  =  2,3,...,M  (2.17) 


The  HOSM,  as  developed  for  wave-wave  interactions,  represents  the  velocity  potential 
at  each  order  as  a  Fourier  series.  The  amplitude  of  each  Fourier  mode  is  determined,  at 
each  order  m,  by  the  Dirichlet  conditions  expressed  above.  In  the  HOSM  developed  for 
solving  wave-body  interaction  problems,  the  velocity  potential  at  each  order  is  expressed  as 
the  sum  of  expansions  of  free  surface  and  body  basis  functions  [23]: 

oo  oo 

$W(x,y,z,t)  =  Y,»(hn)(t)*Fn(x,y,z)  +  Y,°i™Ht)*Bn(x,y,z)  (2.18) 

n=0  n=0 

where  ^in(t)  are  the  Fourier  modal  amplitudes  of  dipoles  distributed  on  the  free  surface  and 
an(t)  are  the  Fourier  modal  amplitudes  of  sources  distributed  on  the  body. 

The  basis  functions,  ^Fn  and  ^Bn,  satisfy  the  Laplace  equation,  spatial  periodicity 
requirements,  and  the  bottom  boundary  conditions.  For  two  dimensional  problems,  for 
example,  the  basis  functions  are  given  as  Fourier  integrals  involving  the  two-dimensional 
periodic  source  potential.  A  derivation  and  explanation  of  the  two-dimensional  periodic 
source  potential  are  available  in  Newman  [34]. 

Boundary  Integral  Equation  Method 

A  boundary  value  problem  governed  by  the  Laplace  equation  can  be  recast  as  a  boundary 
integral  equation  for  the  unknown  velocity  potential  [13].  A  review  of  the  boundary  integral 
equation  method  is  available  in  Hunt  [13]. 

The  boundary  integral  equation  method  (BIEM)  uses  Green's  theorem  to  derive  an 
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integral  equation  for  the  velocity  potential  in  terms  of  distributions  of  singularities  on  the 
boundary  of  a  simply  connected  fluid  domain.  The  velocity  potential  is  expressed  as  a 
function  of  singularity  distributions  on  the  fluid  boundary  by  applying  the  second  form  of 
Green's  Theorem  [12] 


I  f  \JH  -  G- 

J   Jsl    dn         i 


dS  =  0  (2.19) 


where  the  Green  function  G  and  the  velocity  potential  3>  are  solutions  of  the  Laplace 
equation  in  the  fluid  domain. 

If  the  Green  function  is  the  potential  of  a  source  (G  =  logr),  equation  (2.19)  can  be 
replaced  by  [13], 

/  L  m  w  dS<  - 1 L  G™W>dS*  - Tm      (2-20) 

where  T  =  0,  — 27T,  or  —  4n  for  points  outside,  on  the  boundary,  or  inside  the  volume 
defined  by  the  surface  S  =  Sb  U  So  U  Sl  U  Sr  U  Sp  of  Figure  2-1. 

The  choice  of  Green  function  G  is  determined  by  the  boundary  conditions  and  the  type 
of  flow  to  be  modeled.  A  complete  development  of  the  boundary  integral  equation  method 
(BIEM)  is  available  in  Hunt  [13]. 


The  Combined  High-order  Spectral  and  Boundary  Integral  Equation  Method 

The  combined  high-order  spectral  and  boundary  integral  equation  method  (hereinafter  re- 
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ferred  to  as  the  'combined  method')  eliminates  the  limitations  of  the  base  methods  identified 
in  Chapter  1.  In  particular,  arbitrary  body  shapes  are  accommodated  and  all  body  bound- 
ary conditions  are  satisfied  on  the  exact  positions  of  moving  bodies.  The  combined  method 
defines  a  "total  potential",  $r(x,y,z,  £),  as  the  sum  of  a  "spectral  potential",  $sp(x,  y,  z,  t), 
and  a  "body  potential",  $fl(x,t/,  2,  £).  These  potentials  are  defined  by  the  following  rela- 
tionships: 


M 


$T(x,y,z,t)  =  £  *T •    (*.V.*.*)  (2-21) 

m=l 

M 

$sp(x,  y ,  z,  t)  =  £  &W  (x,  y ,  2,  t)  (2.22) 


m=l 

M 


$B(x,y,2,i)  =  £  $T(x,y,z,t)  (2.23) 

m=l 


In  the  combined  method  $t{x,  y,  2,  i)  is  constructed  to  satisfy  the  initial-boundary  value 
problem  specified  above. 

The  function  representing  the  body  potential  is  chosen  to  set  the  body's  contribution  to 
the  potential  on  the  mean  free  surface  to  zero.  For  a  non-lifting  body  the  potential  is  that  of 
a  periodic  array  of  sources  located  on  the  body's  surface  and  its  negative  image.  This  choice 
is  motivated  by  the  sequence  of  boundary  conditions  used  to  determine  the  potential  at  each 
order,  $(m)(x,0, 2,£),  (equations  2.16  and  2.17).  The  sequence  of  boundary  conditions  to 
be  satisfied  on  the  mean  free  surface  is 

4X)  (*,  0,  z,  t)  =  #W  (x,  0,  z,  t)  +  #W  ^  o,  2,  t)  =  $5(x,  0,  z,  t)  (2.24) 
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and 


m_1  ^i  Ft 


771—1        I       C\l 

-£  yr^^^^O^^),        m  =  2,3,...,M         (2.25) 


As  stated  above,  the  body  potential's  contribution  to  the  velocity  potential  on  the  mean 
free  surface  is  zero.  With  this  choice  the  left  hand  sides  of  equations  (2.24)  and  (2.25)  are 
conditions  to  be  satisfied  by  the  spectral  potential.  At  each  order,  the  right  hand  sides 
are  known  Dirichlet  boundary  conditions.  The  coefficients  of  the  Eigenfunction  expansions 
which  represent  the  spectral  potential  are  determined  directly  from  equations  (2.24)  and 
(2.25). 

In  the  combined  method  the  body  boundary  condition  is 


-5—  =  -5^  +  ~^~  =  U  -non  SB{t) 
on  on  on 


(2.26) 


where  U  is  the  body  velocity  and  h  is  the  body  unit  normal  vector. 

In  the  combined  method  —d$sp/dn\sB  is  added  to  the  body's  normal  velocity  to  yield 
the  condition  to  be  satisfied  by  the  singularity  distribution  on  the  body.  For  a  non-lifting 
body  the  condition  is: 


d_ 
dn 


J  js   a(t)G(x,£)dSt 


(m) 


d$sp 


dn 


■\sB+U-h, 


m  =  1 


(2.27) 


d_ 

dn 


J  J    e(OG(x,OdSt 


dn 


\sB  » 


m  =  2,3,  ...,M 


(2.28) 
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where  a(£)  is  the  function  describing  the  source  distribution  on  the  body  surface.  Equa- 
tions (2.27  and  (2.28)  can  derived  from  the  second  form  of  Green's  theorem  [17]  or  by 
directly  assuming  that  the  body  can  be  represented  by  a  distribution  of  sources  along  the 
body  boundary. 

The  coupling  of  the  body  and  and  spectral  potentials  occurs  through  the  body  bound- 
ary condition,  equations  (2.27)  and  (2.28),  and  the  free  surface  boundary  conditions,  equa- 
tions (2.12),  (2.13),  (2.24)  and  (2.25)  .  The  value  of  the  body  potential's  contribution  to 
the  total  potential  on  the  mean  free  surface  is  zero;  the  value  of  the  vertical  component  of 
the  gradient  of  the  body's  potential  on  y  =  0,  d$B/dy\(x,o,z,t)  1S  n°t-  The  body  contributes 
to  the  conditions  to  be  satisfied  on  y  =  0  through  this  quantity. 

Force 

The  force  on  a  body  due  to  the  dynamic  pressure  in  an  ideal,  incompressible,  irrotational 
fluid  is  obtained  by  integrating  the  pressure  over  the  body  surface  [33]: 


'=  /    f    pndS  (2.29) 

J    JsB 

=  f    (    p  (r  x  n)  dS  (2.30) 

J    JSf, 


M 

>sB 


The  useful  forms  of  equations  (2.29)  and  (2.30)  are  derived  by  applying  to  the  unsteady 
form  of  Bernouilli's  equation  either  (1)  the  kinematic  transport  theorem  and  Gauss'  theorem 
or  (2)  the  definition  of  the  exact  differential  of  the  potential.  If  method  (1)  is  applied  the 
resulting  force  expression  is  [33] 

F=AtiiSB*ndS-\(v*v*n)ds+pIL{f^ds    (231) 
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The  following  expression  for  the  force  is  a  result  of  the  application  of  method  (2)  ([16]  and 
[54]) 

F  =  -p^-ff    <t>ndS-l(V<p-V<t>n)dS  +  p  f  f    (U-V0)ndS  (2.32) 

at  J  JsB  2  J  JsB 


The  force  expressions  differ  in  the  third  term  on  their  respective  right  hand  sides.  The 
difference  between  the  expressions  is  explored  in  Chapter  5. 

Equation  of  Motion 

The  equation  of  motion  for  a  rigid  body  is  derived  from  Newton's  Laws.  The  appropriate 
expression  of  dynamic  equilibrium  for  a  direct  simulation  of  the  free  and  forced  motions  of 
a  submerged  body  is  Newton's  Second  Law  [20] 

M  X(i)  =  F(t)  (2.33) 

M  is  the  body  mass  matrix,  X(t)  is  the  body  acceleration  vector,  and  F(t)  is  the  force 
vector.  F(t)  includes  the  wave  and  control  surface  forcing.  The  solution  of  the  equation 
of  motion  for  a  body  enables  simulation  of  free  and  forced  motions.  The  calculation  of 
the  wave  forcing  presently  controls  motion  simulation  computation  speed  and  is  a  principal 
motivation  for  this  thesis. 

A  state  variable  representation  of  the  body  kinematics  makes  direct  use  of  the  solution 
of  equation  (2.33)  to  simulate  body  motions  in  a  time  varying  flow 
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/l+Ol    - 
X{t)dt  (2.34) 

X(t)d<  (2.35) 


Equations  (2.34)  and  (2.35)  are  used  to  update  the  body  velocity  and  position 


X(t  +  6t)  =  X(t)  +  5X(t)  (2.36) 

X(t  +  St)  =  X{t)  +  6X{t)  (2.37) 

In  the  following  section  the  equations  used  to  implement  the  combined  method  in  two 
dimensions  are  presented. 

2.4     Implementation  of  the  Combined  Method  in  Two  Di- 
mensions 

The  combined  high-order  spectral  method  -  boundary  integral  equation  method  formulated 
for  two  and  three  dimensions  has  been  implemented  in  two  dimensions  in  this  thesis.  The 
purpose  of  this  section  is  to  present  the  equations  used  to  carry  out  this  implementation. 

The  chosen  representation  of  the  body  potential  enables  the  initial  condition  for  the 
free  surface  potential,  $s(x,t  =  0),  to  be  used  with  equation  2.24  to  determine  the  initial 
coefficients  of  the  eigenfunctions  which  represent  the  spectral  potential: 


${T](x,0,t  =  0)  =  *£,(ar,0,t  =  0)  +  ^]{x,0,t  =  0)  =  $s(x,t  =  0)  (2.38) 
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With  the  body  potential  represented  as  an  array  of  sources  and  its  negative  image, 

$g)(a:,0,*)=0  (2.39) 

and  the  value  of  the  spectral  potential  on  the  mean  free  surface  is 

$$(x,0,*  =  0)  =  $5(:r,*  =  0)  (2.40) 

The  spectral  potential  eigenfunction  expansion  is  represented  as  a  Fourier  series  and  the 
coefficients  of  each  mode,  <&n3p  (0)  are  determined  by 

N 

*Sj?(*.tM  =  0)  =  £  *gi(0)*n(*,0)  (2.41) 

71=1 

The  representation  of  the  spectral  potential  is  similar  to  that  of  the  HOSM  of  Dommer- 
muth  and  Yue  [8]. 

The  initial  body  velocity,  U(0),  and  the  gradient  of  the  first  order  spectral  potential, 
V<frsp  (x,y),  are  used  with  the  first  order  body  boundary  condition  to  determine  the  first 
order  source  distribution  on  the  body  a^l'(x,y,t): 


d_ 

dn 


f    <jW(x,y,t  =  0)G(x,y,x',y')dSB 
.JsB 


=  C?-n-^g^|SB  (2.42) 
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where  G(x,  y,  x',  y')  is  the  potential  of  a  periodic  array  of  sources  [34]  and  its  negative  image 
reflected  about  y  =  0 


G(WV)    =     ilog(2cosh5(£L^)_2cos^)) 


_     l^2cosh^-2cos^-(2W>)) 


where  h0  is  the  depth  of  the  body's  centroid. 

The  first  order  source  distribution  on  the  body  and  the  first  order  spectral  potential  are 
the  forcing  for  the  solution  of  the  second  order  problem  through  the  boundary  condition 

(2) 

imposed  on  <fr^  (x,0,£): 

42)  (*,  0,  t)  =  *g>  (*,  0,  t)  +  $f  (x,  0,  t)  (2.44) 

As  with  the  first  order  body  potential,  $B' (x,0,t)  =  0.    Therefore  equations  (2.44)  and 

(2) 

equations  (2.25)  are  combined  to  yield  the  condition  to  be  satisfied  by  $sP  (x,0,t): 

^(l,0,t)  =  -,^M-,^lM  (2.45) 

The  second  order  spectral  potential  is  used  in  the  body  boundary  value  problem  to  solve 
for  the  second  order  body  potential. 

-^|5*=-^T|5*  (2-46) 

The  sequence  used  to  determine  the  second  order  source  strengths  and  the  higher  order 
spectral  potentials  follows  that  used  for  the  first  order  quantities  above. 
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With  the  body  and  spectral  potentials  known  to  the  desired  order  M,  the  quantities 
needed  to  compute  the  pressure  and  force  on  the  body  are  known.  For  a  body  free  to 
respond  to  the  waves,  the  force  is  computed  and  used  to  update  the  body's  state  using  the 
methods  described  above. 

The  body  and  spectral  potentials  are  used  with  the  Zakharov  equations  to  update  the 
free  surface  potential  and  free  surface  shape.  The  quantities  needed  are 

dy{x,t)        d$s{x,t)  d$T{x,ri,t) 

— a '      5 '  and     5 v2-47) 

ox  ox  oy 

dr](x,t)/dx  and  d<&s{x,t)/dx  are  calculated  in  spectral  space,  by  the  method  used  by 
Dommermuth,  Yue,  and  Liu  [23].  d$T{x,  V,  t)/dy  is  calculated  by  determining  the  contribu- 
tions of  the  body  and  spectral  potentials  separately.  The  spectral  potential's  contribution 
is  computed  as  in  [23]  and  [8]: 

The  body's  contribution  at  each  order  to  the  fluid  vertical  velocity  on  the  free  surface, 

a££^)  (249) 

can  be  computed  directly  on  the  free  surface  or  by  Taylor  series  expansion  of  the  vertical 
component  of  partial  derivatives  of  the  body  potential  evaluated  on  y  =  0.  For  k  >  1, 
dk$g  /dyk  is  made  computationally  efficient  using  Laplace's  equation  to  reduce  the  number 
of  y  derivatives  required  [23]. 
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Chapter  3 


Numerical  Implementation 


The  initial  boundary  value  problem  defined  in  the  previous  chapter  is  solved  numerically. 
The  purpose  of  this  chapter  is  to  describe  the  numerical  implementation  of  the  mathematical 
formulation.  This  chapter  will: 

•  Present  the  problem  in  the  form  of  systems  of  discretized  equations  to  be  solved  by 
appropriate  numerical  methods; 

•  Analyze  the  accuracy,  consistency,  and  stability  of  the  proposed  numerical  scheme; 
and 

•  Demonstrate  necessary  conditions  for  the  method  to  converge  to  solutions  of  the 
physical  problems  of  interest. 

3.1     System  of  Equations 

The  method  developed  for  solving  the  initial  boundary  value  problem  requires  the  solution 
of  a  system  of  equations  at  each  time  step.  The  system  of  equations  is  determined  by  the 
specified  form  of  the  boundary  conditions  in  the  initial  boundary  value  problem.  The  form 
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of  the  free  surface  boundary  conditions,  the  periodic  nature  of  the  problem  formulation,  and 
the  choice  of  a  boundary  integral  representation  of  the  body  determine  the  characteristics 
of  the  numerical  formulation. 

The  numerical  solution  is  enforced  at  points  along  the  boundary  of  the  computational 
domain  S  =  Sb  U  So  U  Si  U  Sr  U  S>  of  Figure  2-1.  The  solution  is  enforced  at  Nx  evenly 
spaced  points  along  the  horizontal  (x)  axis,  Nb  evenly  spaced  points  around  the  body,  and 
No  evenly  spaced  points  along  the  bottom. 

The  system  of  equations  developed  to  represent  the  boundary  value  problem  at  each 
time  step  is  expressed  symbolically  in  matrix  form  as 


[F]      [LBW]     [Low] 


[P'B]      EM      [Lob] 
[K]      [Lbo]      [Loo] 


*£> 


,(m) 


*(m) 


[R(m)] 


U  ■  n  |s       for    m  =  1 


0  for    m  >  2 


[0] 


(3-1) 


F  is  the  array  of  eigenfunctions  for  the  spectral  potential: 


[F]  = 


eifc0xie|/co|yi  eiA;iiie|fci|yi 


eikox2  e\k0\y2  gifci2;2el^i|y2 


ik„   ii    \k.r    Ivi 

g      Nx     'e'    Nx'"1 


eikNxX2e\kNx\y2 


jk0xN      \k0\yNx        ikixN  tJki\yNt      ...        t*^  xNx    \kNx  \yNx 


(3.2) 
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|f'b  1  represents  the  influence  of  the  spectral  potential  on  the  body: 

[Fb]  =  [V(  )  •  n]5B  (3.3) 


[f'     represents  the  influence  of  the  spectral  potential  on  the  bottom: 

[F'o]  =[*<)•  *\s0  (3-4) 


[LBW]  defines  the  influence  of  the  body  potential  on  the  spectral  potential  on  y=0 

[Low]  defines  the  influence  of  the  bottom  potential  on  the  spectral  potential  on  y=0 

[LBB]  defines  the  body's  self-influence 

[LOB]  defines  the  influence  of  the  bottom  potential  on  the  body 

[LBO]  defines  the  influence  of  the  body  potential  on  the  bottom 

[L00]  defines  the  bottom's  self-influence 

3>sp^mM  is  the  array  of  spectral  potential  modal  amplitudes,  An(t) 

$B^  is  the  body  potential 

$0^m^  is  the  bottom  potential 

Note  that  [LQO]  and  [Low]  need  to  be  computed  only  once. 
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R(m) 


(m). 


is  the  function  representing  the  Dirichlet  boundary  conditions  for  $sp 


m— 1     k    fik 


(3.5) 


where  ^m_fc)  (ar,  0,  £)  =  4™"*°  (*,  0,  t)  +  ^~k)  {x,  0,  t)  +  <^m_fc)  (a;,  0,  t) . 
Note  that  all  the  [L]'  s  contain  the  influences  of  the  sources  and  their  negative  images. 
For  deep  water,  3.1  is  rewritten  as: 


IP]       [Lew] 


[F'a]       [Lbb] 


;(m) 
sp 


* 


(m) 


[r(»>] 


U  ■  n  I ,,       for    m  =  1 

\SB 


0  for    m  >  2 


(3.6) 


The  body  potential  is  represented  as  a  horizontally  periodic  array  of  sources  with  neg- 
ative images  reflected  about  the  mean  free  surface,  y  =  0.  When  present,  the  bottom 
can  be  similarly  represented.  A  consequence  of  this  choice  is  that  the  [LBW]  and  [Low] 
matrices  are  identically  zero.  It  is  noted  that  this  construction  may  appear  similar  to  the 
'pressure  release'  problem,  the  high-frequency  limit  of  the  linear  free  surface  boundary  con- 
dition (cf.  Newman  [33]).  However,  in  that  linearized  problem  the  complete  solution  yields 
$(:r,0,  t)  =  0  whereas  this  is  not  the  case  for  the  problem  at  hand.  The  negative  images 
here  are  used  to  provide  computational  efficiency  without  restricting  the  generality  of  the 
solution.  In  particular  it  permits  fast  calculation  for  the  spectral  potential  at  each  time 
step. 
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3.2     Computational  Steps 

The  time  simulation  of  wave-body  interactions,  developed  here  for  the  deep  water  case  for 
simplicity,  proceeds  as  follows: 

1.  Initial  conditions  (t  =  0)  are  specified  for: 

The  free  surface  profile,  r](x,t); 

Free  surface  potential,  $s(x,t);  and 

Body  position  and  velocity,  SB(x,y,t)  and  \J(x,y,t),. 

2.  The  boundary  value  problem  is  solved  for  the  spectral  potential,  </4p  ,  and  the  source 
strength  distribution  on  the  body,  cr™.  For  the  first  set  of  computations,  at  m  =  1: 

a.  The  spectral  potential  is  determined  by  solving  the  following  for  <ps£  : 


3,{m)     _   \T>{m) 


(3.7) 


b.  With  the  spectral  potential  at  a  given  order  known,  the  source  strengths  are 
determined  by  solving 


[fb]  [•£>]  +  [LBB]  [•!?>]  = 


U  •  n  I  <,       for    m  =  1 


for    m  >  2 


(3.8) 


The  above  steps  (2. a  and  2.b)  are  repeated  up  through  perturbation  order  M.  The 
body  contributes  to  (fyj£   through  its  vertical  derivatives. 
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3.  Compute  the  quantities  necessary  for  the  time  integration  of  the  free  surface  boundary 
conditions.  These  quantities  are: 

$y(X,t](X,t),t)\ 

$f(M);  and 

4.  Integrate  the  free  surface  boundary  conditions,  equations  (2.12)  and  (2.13),  to  update 
the  free  surface  potential,  $s(x,t),  and  free  surface  elevation,  rj(x,t) 

Steps  2,  3,  and  4  are  repeated  at  each  time  step. 

The  following  section  details  the  manner  in  which  each  of  the  above  steps  is  performed. 

3.3     Solution  of  the  Boundary  Value  Problem 

Initial  Conditions-  Free  Surface  Elevation  and  Potential 

The  solution  of  the  nonlinear  wave  body  interaction  problem  requires  the  free  surface 
elevation  and  potential  to  be  accurately  specified.  The  free  surface  profile  and  velocity 
potential  used  here  for  initial  conditions  are  generated  by  Stokes'  expansions.  The  mapped 
equations  of  Schwartz  [41]  are  used  to  specify  a  Stokes'  Wave  of  arbitrary  order  N  [22].  For 
the  typical  wave-body  interaction  analyzed  in  this  thesis  N  =  20  is  sufficient  for  specifying 
the  initial  free  surface  profile  and  velocity  potential.  The  differences  in  the  free  surface 
profile  for  expansions  of  order  2,  10  ,and  20  are  shown  in  Figure  3-1. 
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Free  Surface  Profiles  by  Stokes'  Expansion 
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Figure  3-1:  Initial  free  surface  profile 
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Though  the  calculations  done  here  are  for  Stokes'  waves  as  initial  conditions,  the  equa- 
tions and  solution  methods  apply  to  fully  arbitrary  initial  conditions. 

Influence  Coefficient  Matrices:  [F]  and  [f'bJ 


The  spectral  potentials  at  each  order,  4>%  ,  are  obtained  by  solving  equation  3.7  for  the 
Fourier  modal  amplitudes,  An(t).  In  expanded  form: 


71=0 

$s(x,t)  m=l 

(3.9) 

^  -E^=i1^[^(4r*)(^y,*)+^m~fc)(^^*))]y=0    rn>2 

The  higher  order  spectral  potentials  are  functions  of  the  free  surface  elevation  and 
vertical  derivatives  of  the  lower  order  potentials  evaluated  on  the  mean  free  surface.  The 
vertical  derivatives  are  obtained  by  computing  the  modal  amplitudes  of  the  lower  order 
potentials  by  Fast  Fourier  Transform  (FFT),  multiplying  the  modal  amplitudes  by  |fcn|, 
and  generating  the  derivatives  in  physical  space  by  inverse  FFT: 


dly 

a  n=l 


££_    _    V   Ifr     \lJknXjkn\y 


=  £  \kn\lelknXe\k^\y=0  (3.10) 


The  products  needed  in  the  boundary  conditions  for  the  spectral  potentials,  e.g., 


ri{x)-%-(x,0,t)  (3.11) 

oy 

fvn\ 

are  formed  in  physical  space.  The  operational  count  associated  with  obtaining  each  <psp    is 
0(Nx\ogNx). 
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With  the  Fourier  amplitudes  of  the  spectral  potentials  known,  the  spectral  potential  and 
its  derivatives  can  be  determined  throughout  the  fluid  domain.  In  particular,  the  influence 
of  the  spectral  potential  on  the  body  is  computed  directly  from 


[F'B]  #<»>  =  [v*<™>  •  n 


(3.12) 


A0(t)  ({ikeik°x*)i+  (Nel^)j)  ■  nB 
Ai{t)  ((ikeiklXi)i  +  (Ifeilel*0^)])  ■  nB 

A2{t)  ((ikeik*x*)i  +  (\k2\e^^)S)  ■  nB 

ANx((ikeik»*x>)U(\kNJe^y>)S)-nB 

X{,yi  are  the  locations  of  the  Nb  body  segment  midpoints. 

The  operational  count  associated  with  determining  the  influence  of  the  spectral  potential 
on  the  body  is  0(NxNb)-  The  influence  of  the  spectral  potential  on  the  body  is  combined 
with  the  body  boundary  condition  to  form  the  right  hand  side  for  the  equation  to  be  solved 
for  the  source  distribution  on  the  body: 


[a(xi,yi,t)]  =  [L 


1-1 


BB. 


(3.13) 


(tf-n-^]^)     m  =  l 
"  [P'B]  &  rn  >  2 

Equation  (3.13)  is  solved  by  factoring  the  body-body  influence  coefficient  matrix,  [LBB]~  , 
and  determining  the  source  strengths  by  back  substitution.  The  body-body  influence  co- 
efficient matrix  must  be  factored  only  once,  at  the  first  time  step.  The  operational  count 
associated  with  solving  3.13  at  the  first  time  step  is  0(NB).  The  operational  count  at 
subsequent  time  steps  is  O(NbNb)- 
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Equation  3.13  requires  derivatives  of  the  form 

QkAm-k) 

— ^ (x,(M),    fe  =  l,2,...,m-l  (3.14) 

d<pB  (x,0,t)/dy  is  computed  using  the  analytical  expression  for  the  representation  of  the 
body  as  a  distribution  of  periodic  sources  and  their  negative  images.  Higher  derivatives 
are  found  through  the  use  of  Laplace's  equation  ([8]).  For  example,  [d3(f)B  /dy3)  |  0  is 
required  for  the  solution  of  the  fourth  order  problem.  d<pB  /dy  is  a  harmonic  function, 
therefore 


a2  (o&\  _     a2  (a$ 


(3.15) 


dy2  \    dy   J  dx2  \    dy   J 

d<pB  /dy  is  computed  as  before  and  the  horizontal  derivatives  are  computed  numerically. 
The  symmetry  properties  of  the  body-image  system  force  the  even  numbered  vertical  deriva- 
tives to  be  zero. 

After  the  initial  set  up  effort,  the  total  operational  count  associated  with  solving  the 
boundary  value  problem  to  order  M  at  each  time  step  is 

0(M[Nx  \og(Nx)  +  NXNB  +  NBNB])  (3.16) 

The  total  operational  count  of  the  High  Order  Spectral  Method  (HOSM)  of  Liu,  Dom- 
mermuth,  and  Yue  [23]  is 

0{M[Nx  \og{Nx)  +  NXNB  +  NB  log  NB])  (3.17) 
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As  in  the  HOSM,  typically  Nx  »  Nb-  Thus,  for  large  numbers  of  panels  the  new  method 
is  slower  than  the  full  spectral  method,  but  provides  the  benefit  of  allowing  completely 
arbitrary  body  geometries. 

3.4     Updating  the  Free  Surface  Potential  and  Elevation 

The  Zakharov  equations  (equations  2.12  and  2.13)  are  used  as  time  evolution  equations 
for  the  free  surface  elevation  r)(x,t),  and  the  free  surface  potential  $5,(x,<).  The  time 
integration  requires  (a)  the  accurate  determination  of  the  spatial  derivatives  of  the  free 
surface  elevation  r)x,  free  surface  potential  0f ,  and  velocity  potential  (f>y,  and  (b)  a  stable 
and  accurate  time  integration  scheme. 

Horizontal  Derivatives  of  r](x,t)  and  4>s(x,t) 

The  free  surface  elevation  and  free  surface  potential  are  defined  at  evenly  spaced  points 
along  the  x  axis.  The  horizontal  derivatives  of  these  quantities  are  most  efficiently  performed 
in  Fourier  space: 


^jM     =     £  Bn{t)ikne*«x  (3.18) 

n=0 

^1^-     =     £  A*(i).-*„e"»*  (3.19) 

n=0 


where  Bn(t)  and  A% (t)  are  the  Fourier  amplitudes  of  the  free  surface  elevation  and  potential 
respectively.  The  operational  count  associated  with  computing  these  spatial  derivatives  is 
0(NX  \ogNx).  Products  involving  these  quantities  are  computed  in  physical  space. 

The  Vertical  Velocity  of  the  Fluid  on  the  Free  Surface 
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The  vertical  velocity  of  the  fluid  on  the  free  surface  is  the  sum  of  the  contribution  from 
the  spectral  potential  and  the  body  potential.  The  contribution  from  the  spectral  potential 
is  determined  by  computing  the  derivatives  in  Fourier  space  [8]: 


dk+i 


\k+i 


O^Ti*n(x,y)\y=o  =  \knrl*n(x,y)\y=0  (3.20) 


The  Fourier  amplitudes  are  found  by  FFT  and  the  appropriate  power  of  the  wave  number 
is  multiplied  by  the  series  representation  of  the  function,  and  the  required  derivative  of  the 
function  is  computed  by  the  inverse  FFT.  The  operational  count  associated  with  computing 
these  derivatives  is  0(NX  log  Nx). 

The  contribution,  at  each  order,  of  the  body  potential  to  the  fluid  vertical  velocity  on 
the  free  surface  is  computed  using  the  derivative  of  analytical  expression  for  the  potential 
due  to  the  body  and  its  negative  image. 

The  operational  count  associated  with  computing  the  vertical  derivatives  is  0(NxNb)- 
The  total  operational  count  associated  with  computing  the  terms  required  for  solution  of 
the  kinematic  and  dynamic  free  surface  boundary  conditions  is  0{M[Nx  log  Nx  +  NxNb  + 
NBNB}). 

The  above  steps  provide  the  information  required  to  compute  dr](x,  t)/dt  and  d(f)S(x,  t)/dt 
at  each  time  step.  In  the  next  section  the  method  used  to  integrate  the  free  surface  evolution 
equations  in  time  is  presented. 

Time  Integration 

The  integration  scheme  used  to  evolve  the  free  surface  shape  and  free  surface  potential 
is  a  fourth  order  Runge-Kutta  (RK4)  scheme.  The  stability  properties  of  this  method  as 
applied  to  the  non-linear  wave  body  interaction  problems  were  developed  by  Dommermuth 
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[7]  and  used  by  Liu  [22]  in  his  extension  of  Dommermuth's  work.  The  important  points 
from  Dommermuth  are: 

1.  RK4  is  conditionally  stable; 

2.  The  Courant  condition  for  RK4,  based  on  von  Neumann  stability  analysis  of  the 
linearized  free  surface  conditions,  is  g  kgrid&t2  <  8,       kgrid  =  2n/Ax;  and 

3.  RK4  is  mildly  dissipative. 

Analyses  of  other  time  integration  schemes  for  numerical  simulations  of  wave-wave  and 
wave-body  interactions  are  found  in  Kring  [20]  and  Dommermuth  [7]. 

Considerable  information  is  available  concerning  the  development  and  growth  of  "saw- 
tooth instabilities"  in  time  simulations  of  free  surface  flows.  Olmez  [37]  and  Dommermuth 
[7]  provide  comprehensive  background  discussions  of  the  saw-tooth  instability.  Longuet- 
Higgins  and  Cokelet  [25]  were  the  first  to  report  the  presence  of  the  sawtooth  instability 
and  contended  that  its  presence  was  partly  physical.  Other  investigators  reporting  sawtooth 
type  instabilities  are  Baker  et  al.  [1],  Lin  et  al.  [21],  and  Sen  et  al.  [42]. 

Dommermuth  [7]  and  Roberts  [39]  stated  that  the  sawtooth  instabilities  are  numerical. 
Olmez  [37]  stated  that  the  cause  of  the  instability  is  unknown. 

The  discretization  of  the  physically  continuous  problem  introduces  a  sawtooth  wave 
form  at  the  highest  resolvable  wave  number,  k  =  27r/Ax  (Figure  3-2). 

The  coupling  provided  by  the  method  of  solution  outlined  in  the  previous  chapter  ensures 
that  all  modes  at  all  orders  are  coupled.  Instabilities  at  any  mode  are  coupled  to  all  other 
modes.  The  presence  of  physical  or  numerical  stabilizing  mechanisms  determines  whether  a 
given  frequency  is  stable.  In  the  chosen  mathematical  formulation  of  the  physical  problem 
the  fluid  is  assumed  to  be  ideal.  Therefore  no  physical  dissipative  mechanism  exists.  The 
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The  value  of  the  eigenfunction  at  the  highest  wave  number 

e        ,  k=n/Ax 


II  \         Ax 


V 


2Ax  v        3Ax 


Figure  3-2:  Sawtooth  at  Highest  Wavenumber 

numerical  scheme  must  provide  the  dissipation  absent  in  the  chosen  physical  model.  RK4 
is  mildly  dissipative. 

The  level  of  accuracy  of  the  spatial  derivatives  evaluated  on  the  free  surface  influences 
the  onset  of  sawtooth  instabilities  (Roberts  [39]  and  Dommermuth  [7]).  Derivatives  of  the 
potential  of  each  order  m  are  of  the  form 


dT 


dy> 


■*n  =  i*»r*» 


(3.21) 


where  n  is  the  Fourier  mode.  The  magnitude  of  the  derivative  is  proportional  to  (kn)m.  The 
error  at  a  given  mode  and  order  is  amplified  by  this  factor  [7].  If  dissipation  is  not  provided 
the  growth  of  this  error  can  be  unbounded.  Longuet-Higgins  [25],  Dommermuth  [7],  and  Liu 
[22]  all  use  a  form  of  smoothing  to  further  control  high  wave  number  instabilities.  In  this 
thesis  an  ideal  low  pass  filter  is  used  to  remove  the  highest  wave  number  modes.  The  cutoff 
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wave  number  is  constrained  by  physical  considerations.  The  filter  must  allow  physically 
relevant  waves  to  develop,  Kutoff/Dbody  <^  1  and  must  not  remove  appreciable  energy  from 
the  system.  The  filter  is  implemented  by  transforming  the  free  surface  elevation  and  free 
surface  potential  into  the  frequency  domain  by  the  FFT.  All  modes  above  the  cutoff  wave 
number  are  removed,  and  the  filtered  free  surface  elevation  and  potential  are  generated  by 
inverse  FFT. 

3.5     Convergence 

The  Lax-Richtmyer  equivalence  theorem  states  that  a  consistent  approximation  to  a  well- 
posed  linear  problem  is  stable  if  and  only  if  it  is  convergent  [10].  First,  it  is  necessary  to 
demonstrate  that  the  numerical  solution  method  is  consistent.  The  consistency  is  demon- 
strated by  conservation  of  volume  and  conservation  of  energy.  Confidence  in  converged 
results  is  provided  by  comparison  of  converged  results  with  experimental  results  and  results 
generated  by  other  numerical  methods  (Cointe  [6])  . 
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Energy  and  Volume  Conservation 

Volume  conservation  requires 


tSF{i 
Energy  conservation  requires 


/        r](x,  t)  dx  =  Constant  (3.22) 

JSr(t) 


"IIIv{lv2  +  s«)dV-C  (323) 

A  convenient  form  of  the  total  energy  can  obtained  by  making  use  of  the  kinematic  infor- 
mation, U  =  V0  and  V20  =  0  and  the  first  form  of  Green's  theorem 

I  [<t>V2<f>  +  (V</>)2]  dv  =  /  n  ■  <j)V(j>ds  (3.24) 

The  above  information  can  be  used  to  express  the  kinetic  energy  T  as  (Milder  [28]  and 
Kinsman  [18] 

T=-pfv<t>-  V<f>dv  =  7,p[  ^V<£  '  ds  ~  2P  J  ^  dx  ^3-25^ 

The  expression  used  for  computing  the  total  energy,  £,  in  the  system  is 

C  =  -pj<f>SVt  dx+-pg   I V  dx  (3.26) 

The  total  energy  is  computed  at  each  time  step. 
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Convergence-  Discussion 

The  convergence  properties  of  the  combined  HOSM-BIEM  method  are  functions  of  the 
following  sources  of  error  [7]  : 

1.  Finite  numbers  of  body  segments  Nb,  Fourier  modes  Nx,  and  perturbation  order  M  ; 

2.  Time  integration  scheme; 

3.  Finite  computational  domain; 

4.  Filter  cutoff  wave  number  Nfuter\  and 

5.  Aliasing. 

The  first  four  sources  of  error  require  no  explanation.  The  convergence  properties  associated 
with  each  will  be  examined  in  the  subsequent  section. 

Solutions  associated  with  wave  numbers  above  the  highest  grid  wave  number  (k  = 
27r/A:r)  can  not  be  distinguished  from  waves  of  lower  frequency  [46].  This  is  aliasing  and 
it  occurs  is  when  products  are  formed  in  the  spectral  method,  e.g.,  tjx t}x 


Nx  N; 

E  Vn(t)eihnXikn  I 
n-0  n=0 


Y^Mt^^flMt^^^-^M2^^  +■■■+-*£.  {vNx)2e2kN*     (3-27) 


Products  are  made  alias  free  using  the  method  developed  by  Orszag  [38]. 
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Convergence-  Demonstration 

The  convergence  properties  of  the  HOSM-BIEM  are  demonstrated  in  this  section. 

Nomenclature: 
M  =  Perturbation  order 
Nb  =  Number  of  body  segments 
Fy  =  Steady  vertical  force 
p  =  Fluid  density 
g  =  Gravitational  constant 
A  =  Wave  amplitude 
k  =  Wave  number 
R  =  Body  radius 
H  =  Distance  from  mean  free  surface  to  center  of  body 

Nw  =  Domain  size;  the  number  of  waves  of  the  fundamental  wave  length  in  the 
domain 

Nx  =  Number  of  Fourier  modes;  equal  to  the  number  of  segments  along  the  x  axis 

T  =  Period  of  the  fundamental  wave 

Ts  =  Simulation  time 

A£  =  Time  step  size 

N futer  =  Fourier  mode  for  low  pass  filter  cut-off 
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NB 

M  =  2 

M  =  3 

M  =  4 

64 

0.2807 

0.2799 

0.2798 

128 

0.2783 

0.2775 

0.2775 

256 

0.2768 

0.2760 

0.2760 

Table  3.1:  Convergence  of  the  vertical  drift  force,  Fy/pgA2,  with  body  segments  Nb  and 
order  M.  kA  =  0.04,  kR  =  0.4,  and  H/R  =  2;  and  Nw  =  16,  Nx  =  64NW,  NfMer  =  24NW 
for  NB  =  64  and  NB  =  128,  NIilter  =  32ATW  for  NB  =  256,  T/A*  =  128,  Ts  =  5T. 


Table  3.1  shows  the  convergence  of  the  vertical  drift  force  with  increasing  number  of 
body  segments  Nb  and  perturbation  order  M.  The  convergence  rate  appears  to  be  better 
than  linear  with  increasing  number  of  body  segments,  and  exponential  with  increasing 
perturbation  order.  The  full  spectral  method  of  Dommermuth  converges  exponentially 
with  increasing  number  of  body  modes  (segments).  The  combined  method  of  this  thesis 
uses  a  boundary  integral  representation  for  the  body's  contribution  to  the  velocity  potential. 
Therefore  the  convergence  rate  is  not  expected  to  match  that  of  the  full  spectral  method. 
The  decreased  convergence  rate  is  a  penalty  associated  with  the  improved  flexibility  of  the 
present  method. 

Figure  3-3  shows  the  behavior  of  the  solution  for  the  vertical  force  as  the  value  of  the 
low-pass  filter  cut-off  wave  number  is  varied.  The  cut-off  wave  number  corresponds  to  waves 
whose  lengths  vary  from  2/3  to  1/48  of  the  wavelength  of  the  fundamental  of  the  initial 
Stokes  wave.  For  cut-off  wave  numbers  corresponding  to  waves  shorter  than  1/4  of  the 
fundamental,  there  is  no  appreciable  change  in  the  vertical  drift  force. 

Table  3.2  shows  the  convergence  of  the  vertical  drift  force  with  increasing  number  of 
Fourier  modes  and  perturbation  order.  The  convergence  rate  appears  to  be  exponential  with 
increasing  number  of  Fourier  modes  and  exponential  with  increasing  perturbation  order. 
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Figure  3-3:  Behavior  of  the  vertical  drift  force,  Fy/pgA2,  with  varying  cut-off  wave  number 
Nfuter-  kA  =  0.04,  kR  =  0.4,  and  H/R  =  2;  and  NB  =  256,  Nw  =  16,  Nx  =  64NW, 
T/At  =  128,  M  =  4,  Ts  =  ST. 


Table  3.3  shows  how  the  convergence  of  the  vertical  drift  force  varies  with  increasing 
wave  slope.  The  convergence  rate  with  increasing  number  of  Fourier  modes  decreases  as 
the  wave  slope  increases.  Nx  =  32iV\v  Fourier  modes  was  insufficient  for  kA  =  0.12.  The 
number  of  Fourier  modes  was  insufficient  to  resolve  the  local  steep  waves  in  the  vicinity  of 
the  body  and  the  numerical  solution  broke  down  after  2  wave  periods. 

Table  3.4  shows  the  convergence  of  the  method  with  respect  to  time  step  size.  The 
behavior  of  the  solution  with  time  step  size  is  consistent  with  the  findings  of  Korsmeyer 
[19].  For  solutions  that  are  converged  with  respect  to  the  spatial  variables,  there  is  a  wide 
range  of  solutions  for  which  changes  in  time  step  size  do  not  affect  the  accuracy  of  the 
solution.  For  larger  time  step  sizes  the  solution  error  grows  rapidly  and  is  unstable  [19]. 
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e 

Nx 

M  =  2 

M  =  3 

M  =  4 

0.04 

32  Nw 

0.2769 

0.2758 

0.2760 

64  Nw 

0.2768 

0.2760 

0.2760 

128  Nw 

0.2768 

0.2760 

0.2760 

Table  3.2:  Convergence  of  the  vertical  drift  force,  Fy/pgA2,  with  Fourier  modes  Nx  and 
order  M.  kR  =  0.4,  and  H/R  =  2;  and  Nw  =  16,  Nfilter  =  24NW,  NB  =  256,  T/At  =  128, 
Ts  =  5T. 


€ 

32  Nw 

64  Nw 

128  Nw 

0.04 
0.08 
0.12 

0.2760 

0.2760 

0.2760 

0.2760 

0.2706 

0.2706 

- 

0.2692 

0.2640 

Table  3.3:  Convergence  of  the  vertical  drift  force,  Fy/pgA2,  with  Fourier  modes  Nx  as  wave 
slope  kA  increases.  M  =  4,  kR  =  0.4,  and  H/R  =  2;  and  Nw  =  16,  Njuter  =  24Nw, 
NB  =  192,  T/At  =  128,  Ts  =  ST. 


Figure  3-4  shows  the  behavior  of  the  HOSM-BIEM  method  as  the  domain  size  is  varied. 
The  purpose  of  evaluating  the  behavior  of  the  solution  as  the  domain  size  varies  is  to 
determine  the  influence  of  the  periodic  problem  formulation  on  the  ability  to  achieve  and 
sustain  limit-state  behavior  before  "reflections"  from  the  periodic  boundaries  effect  the 
solution.  Figure  3-4  shows  that  limit  state  behavior  is  reached  within  3  periods,  Ts  =  3T, 
of  simulation  time.  The  effect  of  "reflections"  for  Nw  =  8  is  evident  in  Figure  3-4. 
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T/At 

M  =  2 

M  =  3 

M  =  4 

32 

0.27826 

0.2775 

0.2774 

64 

0.2783 

0.2775 

0.2775 

128 

0.2783 

0.2775 

0.2775 

Table  3.4:  Convergence  of  the  vertical  drift  force,  Fy/pgA2,  with  time  step  size  T/At  and 
order  M.  kA  =  0.04,  kR  =  0.4,  and  H/R  =  2;  and  Nw  =  16,  NfMer  =  2WW,  NF  =  G4NW, 
NB  =  128,  Ts  =  5T. 


0.290 


0.265 


Figure  3-4:  Behavior  of  the  vertical  drift  force,  Fy/pgA2,  with  varying  domain  size  Nw 
where  Nw  is  the  ratio  of  the  domain  size  to  the  fundamental  wavelength.  kA  =  0.04, 
kR  =  0.4,  and  H/R  =  2;  and  NB  =  128,  Nx  =  64iVV,  Nfilter  =  24NW,  T/At  =  128, 
M  =  4 
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Chapter  4 


Wave  Diffraction  by  Submerged 
Cylinders 


The  purpose  of  this  chapter  is  to  present  the  results  for  the  analysis  of  the  diffraction  of 
incident  waves  by  submerged  circular  cylinders.  The  diffraction  of  waves  by  cylinders  of 
circular  and  'pontoon'  cross  sections  are  studied.  The  pontoon  case  is  analyzed  to  demon- 
strate the  applicability  of  the  method  to  arbitrary  geometric  shapes.  Figure  4-1  defines  the 
problem  parameters. 

The  oscillatory  and  steady  forces  for  a  range  of  problem  parameters  are  presented  and 
compared  with  related  results  generated  from  the  other  solution  methods  and  experimental 
results  discussed  in  Chapter  1. 

4.1     Diffraction  by  a  Cylinder  of  Circular  Cross  Section 

The  forces  resulting  from  the  interaction  of  a  fully  submerged  circular  cylinder  with  Stokes' 
waves  are  analyzed  in  this  section.  The  forces  on  the  cylinder  are  evaluated  as  functions  of: 
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k  =  27i/X 


H 


\* 2R 


2R 


Circular  cross  section 


Pontoon  cross  section 


Figure  4-1:  Diffraction  problem  parameters 

•  Keulegan-Carpenter  number,  Kq  =  ire~kHA/R  ; 

•  Wave  slope,  kA  ;  and 

•  Submergence  ratio,  h  =  H/R. 

Although  a  "natural"  nondimensional  variable  is  kR,  the  variable  Kc,  which  is  equiv- 
alent in  the  presence  of  the  other  nondimensional  variables,  is  used  instead  to  make  the 
results  directly  comparable  to  others  in  the  literature. 

The  forces  analyzed  are: 

•  First,  second,  and  third  harmonics  of  the  vertical  and  horizontal  forces,  and 

•  Steady  vertical  and  horizontal  forces. 
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Method  of  Computing  Forces 

Time  simulations  of  the  initial  boundary  value  problems  are  performed  long  enough 
for  the  solution  to  achieve  limit-state  behavior  and  terminated  when  "reflections"  from  the 
domain  boundaries  begin  to  influence  the  solution.  Limit  state  behavior  is  typically  reached 
within  two  periods.  The  length  of  time  the  limit-state  behavior  is  sustained  is  determined 
by  domain  width,  Nw  (see  figure  3-4).  For  the  diffraction  problem  the  forces  are  computed 
by  applying  Equations  2.29  and  2.30  to  the  time  history  of  the  sum  of  the  potentials  at 
each  order,  fy^  > m  =  1?  2, ... ,  M,  to  generate  a  force  time  history.  For  stationary  bodies 
the  equations  reduce  to 


F(t)  =  -pf    p(t)ndS  =  -pjs     [5^3  +  lv<pT{t)  ■  V<pT(t) 


ndS 


(4.1) 


where  (f>r{t)  is  the  total  potential  on  the  body's  surface. 
For  the  circle  Equation  4.1  becomes 


F(t)  =-pf 
JSB 


d4>T(t)     l/i  d<pT(t) 


')' 


dt 


-(- 
2  \r0 


dd 


ndS 


(4.2) 


where  r0  is  the  circle  radius. 
Numerically, 


NB 

F{t)  =  -pY, 


i=l   L 


d<f>rT(t) 

dt 


+ 


-(- 
2\r0 


1  d<f>iT(t) 


89 


nidS 


(4.3) 


d(piT(t)/dt  and  (d(f>T(t)/d9)i  are  computed  by  temporal  and  spatial  central  differencing, 
respectively.  The  potential  and  pressure  along  each  segment  are  assumed  to  be  constant. 
Therefore  the  integrated  quantity  at  each  time  step  is  computed  as  the  simple  sum, 
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NB 

^PifiidSi  (4.4) 

i=l 
The  force  at  each  order  m  is  related  to  the  nth  harmonic  coefficient  of  the  force: 


/   x       l    rt+T 
Fn  =  F(m)~-/        F{t)e?m{jjtdt-  n  =  m  =  l,2,...,M  (4.5) 

T  Jt 

with  j  =  \/— I. 

The  forces  computed  by  the  combined  method  are  compared  with  solutions  generated 
by  the  other  methods  described  in  Chapter  1.  As  noted  by  Liu,  et  al.  [23],  there  is  no  direct 
relationship  between  the  F^  terms  of  the  nonlinear  solution  and  the  Fn  =  £m=i  em4% 
terms  of  a  perturbation  series  solution  generated  by  a  boundary  condition  perturbation 
frequency-domain  approach.  Liu  et  al.  demonstrate  that  relationships  between  the  nth- 
harmonic  of  the  nonlinear  solution,  <j>nM>  =  Em=i  4>n  •>  n  =  1>  2,  •••,  Nx  are  related  to  the 
frequency  domain  amplitudes,  &n>  by 


4>W  =  0<2>  +  0(€3)  (4fi) 

0(n)=0<M>+o(en+l)  (4J) 

Numerical  Results 

In  Figure  4-2  the  first  harmonic  of  the  horizontal  diffraction  force  on  the  submerged 
circular  cylinder  is  shown  as  a  function  of  the  Keulegan-Carpenter  number.  The  results  are 
presented  together  with: 

•  The  experiments  of  Chaplin  [3]  and  Miyata  [30]; 
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The  linear  theory  of  Ogilvie  [35]; 


The  High-order  Spectral  Results  (HOSM)  of  Liu,  Dommermuth,  and  Yue  [23];  and 


The  Mixed  Eulerian-Lagrangian  (MEL)  method  of  Cointe  [6]. 
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Figure  4-2:  The  first-harmonic  of  the  horizontal  diffraction  force.  Experiment  (x);  linear 
result(— );  HOSM(D);  MEL(A);  present  method  (•).  (kR  =  0.21,  H/R  =  2.0.) 


The  agreement  between  the  method  developed  in  this  thesis  and  the  other  numerical 
methods  for  this  case  helps  to  confirm  the  theoretical  and  numerical  validity  of  the  new 
method. 

In  Figure  4-2  the  experimental  results  differ  appreciably  from  the  theoretical  and  nu- 
merical results  at  values  of  the  Keulegan-Carpenter  number  greater  than  0.5.  The  difference 
is  due  to  flow  separation,  a  fluid  behavior  not  accounted  for  in  the  potential  flow  equations 
used  in  each  of  the  theoretical  and  numerical  methods  cited. 
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The  lack  of  agreement  between  the  experimentally  generated  results  and  the  results 
generated  by  the  numerical  methods  is  due  to  lift  oscillating  at  the  fundamental  frequency 
in  antiphase  with  the  inertia  force  (Chaplin  [3]).  Chaplin's  results  were  consistent  with  the 
analytical  results  of  Longuet-Higgins  [24]  who  computed  the  magnitude  of  the  circulation 
associated  with  the  oscillatory  flow  around  the  cylinder.  The  decrease  in  vertical  force 
with  increasing  Keulegan-Carpenter  number  is  due  to  the  (negative)lift  asscociated  with 
the  circulation. 

The  zeroeth,  second,  and  third  harmonics  of  the  horizontal  force  are  shown  in  Fig- 
ures 4-3,  4-4,  and  4-5  along  with  the  results  generated  by  experiments,  second-order  theory 
(Ogilvie  [35]),  perturbation  theory  (Vada  [48]),  and  HOSM  and  MEL  computations.  The 
results  of  the  theoretical  and  numerical  methods  agree  with  the  experimental  results  over 
a  wide  range  of  Keulegan-Carpenter  numbers. 
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Figure  4-3:  The  steady  vertical  diffraction  force.   Experiment  (x);  MEL(O);  second-order 
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The  drift  forces  as  functions  of  wave  slope  are  shown  in  Figure  4-6.  The  results  of  the 
present  method  are  compared  with  those  generated  by  the  HOSM  for  the  vertical  force 
and  the  HOSM  and  MEL  for  the  horizontal  force.  The  normalized  vertical  drift  force  is 
independent  of  the  wave  slope.  The  normalized  horizontal  drift  force  is  shown  to  decrease 
in  magnitude  as  wave  slope  increases.  The  negative  horizontal  drift  force  was  studied 
extensively  by  Dommermuth  [7]  and  Liu  et  al.  [23]  and  is  not  a  principal  focus  of  this 
thesis.  Cointe  [6]  also  computed  a  negative  horizontal  drift  force.  In  Cointe's  work  the 
negative  horizontal  drift  force  is  only  present  for  wave  steepnesses  close  to  the  limiting 
steepness  of  the  method. 

As  wave  slope  becomes  vanishingly  small,  the  fourth  order  force  should  approach  a 
constant.  As  seen  in  Figure  4-6  none  of  the  numerical  methods  capture  this  result.  The 
errors,  however,  for  the  high  order  spectral  method  and  the  present  method  are  small. 
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Figure  4-6:  Drift  forces  as  functions  of  wave  slope,  kA.   Fy  is  the  vertical  force;  Fx  is  the 
horizontal  force.  HOSM(D);  MEL(A);  present  method(x).  (kR  =  0.4,  H/R  =  2.0.) 
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The  dependence  of  the  drift  forces  on  submergence,  H/R,  is  shown  in  Figures  4-7  and  4- 
8.  The  experimental  results  are  from  Miyata  [30].  As  the  body  approaches  the  free  surface 
the  effects  of  wave  breaking  become  important  and  the  present  method  can  not  account 
for  overturning  waves.  This  limitation  also  exists  with  the  MEL  method  of  Cointe  and  the 
HOSM  of  Liu  et  al. 
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4.2     Diffraction  by  a  Cylinder  of  Pontoon  Cross  Section 

The  ability  of  the  combined  method  to  compute  forces  on  bodies  of  arbitrary  geometry 
is  demonstrated  by  computing  the  forces  resulting  from  the  interaction  of  a  cylinder  of 
pontoon  cross  section  with  incident  waves.  The  pontoon  shape  is  shown  in  Figure  4-1.  The 
cylinder's  geometric  definition  is  the  one  used  by  Vada  [48]: 


cos  9  —  a  cos  30 


x  = 


y  =  b 


l-a 
sin  6  +  a  sin  30 


l-a 


(4.8) 


For  the  shape  evaluated  in  this  thesis,  a  =  0.1  and  b  =  0.5. 


Method  of  Computing  Forces 


For  the  pontoon  Equation  4.1  becomes 


F{t)  =  -pf 

JSg 


d<f>T(t)     l  (d<f>T{t) 

dt         2  V     ds 


ndS 


Numerically, 


NB 

F(t)=-PY: 


i=i 


d<f>iT(t)  i  (d<f>iT(t) 

dt  2  V      ds 


n;  dSi 


(4.9) 


(4.10) 


d(piT(t)/dt  and  (d<f)T(t)/ds)i  are  computed  by  temporal  and  spatial  central  differencing, 
respectively. 

The  harmonics  of  the  force  time  histories  are  determined  using  the  method  presented  in 
the  previous  section. 
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Numerical  Results 

Time  simulations  of  the  interaction  of  the  pontoon  with  waves  were  performed  and  the 
resultant  forces  compared  with  those  of  Vada[47].  The  zeroeth,  first,  and  second  harmonics 
of  the  horizontal  and  vertical  forces  as  functions  of  the  body  submergence,  h  =  H/R,  and 
nondimensional  wave  number,  K  =  R/(u)2g)  are  presented  in  Figures  4-9,  4-10,  and  4-11. 
The  results  for  the  two  methods  are  in  agreement  for  the  range  of  depths  and  wave  numbers 
evaluated. 


Figure  4-9:  First-  and  second-order  horizontal  oscillatory  forces  on  the  pontoon,  (a)  H/R  = 
1.0;  (b)  H/R  =  1.25;  (c)  H/R  =  1.5;  Vada  (----,  \Fxl/pgRA\;  -,  \Fx2/pgA2\);  present 
method  (n,\Fxl/pgRA\  ;  O,  \Fx2/pgA2\) 
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1.0 


Figure  4-10:  First-  and  second-order  vertical  oscillatory  forces  on  the  pontoon,  (a)  H/R  = 

1.0;  (b)  H/R  =  1.25;  (c)  H/R  =  1.5;  Vada  ( ,  \Fyl/pgRA\;  -,  \Fy2/pgA2\);  present 

method  (D,|Fyl/p^^|  ;  O,  \Fy2/pgA2\) 
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h=  1.0 


Figure  4-11:  Steady  vertical  force  on  the  pontoon,    (a)  H/R  =  1.0;  (b)  H/R  =  1.25;  (c) 
H/R  =  1.5;  Vada  (-,  \Fy/pgA2\)\  present  method  (D,  \Fyi/pgA2\) 
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4.3     Summary 

The  ability  of  the  combined  high-order  spectral-  boundary  integral  equation  method  to 
estimate  the  forces  due  to  the  interaction  of  a  stationary  two  dimensional  cylinder  of  arbi- 
trary cross  section  with  waves  has  been  demonstrated.  The  method  enables  second  order 
and  higher  forces  to  be  computed  with  accuracy  consistent  with  the  accuracy  of  high-order 
spectral  and  Mixed  Eulerian-Lagrangian  methods  for  similar  problem  parameters.  The 
present  method  computes  free  surface  quantities  faster  than  an  MEL  and  can  accommodate 
more  general  body  shapes  than  the  high-order  spectral  method. 

In  the  next  chapter,  the  combined  method  will  be  applied  to  the  radiation  problem  for 
a  fully  submerged  cylinder  undergoing  large  amplitude  forced  oscillatory  motion. 
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Chapter  5 


Wave  Radiation  by  a  Circular 
Cylinder  Undergoing  Large 


Amplitude  Oscillatory  Motion 


The  purpose  of  this  chapter  is  to  demonstrate  the  application  of  the  combined  high  order 
spectral-boundary  integral  equation  method  to  the  problem  of  a  fully  submerged  circular 
cylinder  undergoing  large  amplitude  forced  oscillatory  motion  beneath  an  initially  quiescent 
free  surface.  Figure  5-1  defines  the  problem  parameters. 

The  cylinder  is  forced  to  oscillate  in  heave  and  in  circular  orbital  motion.  The  free 
surface  profiles  and  the  oscillatory  and  steady  forces  resulting  from  the  forced  motions  are 
presented  and  compared  with  results  generated  by  the  other  solution  methods  described  in 
Chapter  1. 
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Figure  5-1:  Radiation  Problem  Parameters 

5.1     Radiation  Problem  for  a  Circular  Cylinder  in  Heave 

The  forces  and  free  surface  profiles  resulting  from  the  interaction  of  a  fully  submerged 
cylinder  forced  to  oscillate  in  heave  with  an  initially  quiescent  free  surface  are  presented  in 
this  section.  The  forces  are  evaluated  as  functions  of: 

•  The  nondimensional  oscillation  amplitude  ,  £/R; 

•  The  submergence  ratio,  h  =  H/R;  and, 

•  The  body  radius  to  wavelength  ratio  (x  2n),  kR  =  2'kR/X. 

The  forces  analyzed  are: 

•  The  first  and  second  harmonics  of  the  vertical  force;  and, 

•  The  steady  vertical  force. 
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Method  of  Computing  Forces 

The  mathematical  formulation  developed  in  this  thesis  solves  the  initial  boundary  value 
problem  in  an  earth  fixed  inertial  reference  system.  For  the  radiation  problem  the  expression 
for  the  force  is  developed  from  Equation  4.1: 


nt)  =  -rJlsm  +  \v*r-V*r 


ndS  (5.1) 


isB 

Where  <f>T  is  the  total  potential  on  the  body  surface. 

The  chosen  numerical  formulation  provides  information  to  compute  d<pr/dt\s  directly. 
However,  the  body  is  moving  so  that  dcfr/df  ^  dfyrjdt.  The  exact  differential  of  the 
potential  is  used  to  develop  an  expression  for  d(pr/dt: 

d(f>r  _  d(pr      dx  d^r      dy  dcpr  _  d$r       -  ,      . 

dt    ~    dt        dt  dx        dt  dy    ~    dt  ^  [      ' 

Equation  (5.2)  is  combined  with  Equation  (5.1)  to  provide  the  expression  used  for  com- 
puting the  force  on  the  moving  body  [16]: 

F  =  -pjjs    [^  -  i  (V0r  •  V0t)  +  (p  ■  V0t)     ndS  (5.3) 

In  Chapter  2,  two  possible  formulations  for  computing  the  force  on  a  moving  body,  equa- 
tions (2.31)  and  (2.32).  Wu  [54]  performed  force  calculations  for  a  surface-piercing  vessel 
using  both  equations  and  concluded  that  the  two  equations  produce  different  answers.  The 
correct  expression  for  computations  involving  bounded  fluids  is  equation  (2.32).  Equation 
(2.31)  is  appropriate  if  the  body  is  surround  by  a  fixed  control  surface  on  which  U  •  n  =  0 
[33]. 
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The  time  domain  simulation  provides  the  force  time  history  in  the  same  fashion  as  in 
the  diffraction  problem  presented  earlier  in  this  thesis.  The  nth  harmonic  coefficients  of  the 
force  are  computed  and  related  to  the  force  at  each  order  m  using  the  method  presented  in 
Chapter  4. 
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Numerical  Results 

The  results  for  the  radiation  problem  are  compared  with: 

•  The  body-exact  linearized  free  surface  formulation  of  Wu  [52]  for  large  amplitude 
oscillations; 

•  The  nonlinear  BIEM  of  Silva  and  Peregrine  [43]; 

•  The  linear  analytic  solution  of  Ogilvie  (1963)  for  first  order  forces  and  the  steady 
component  of  the  second  order  force. 

Trajectory 

The  cylinder  is  started  from  rest  and  the  amplitude  of  the  heave  oscillation,  £(£),  is 
increased  from  zero  to  its  limit  state  value  using  the  following  function  [26]: 

a*)=(l-exp-Qt)£maisinM)  (M 

where  u  is  the  frequency  of  oscillation  and  a  is  a  startup  constant.  Equation  (5.4)  is  used 
to  provide  a  smooth  startup. 

Radiation  Condition-  Numerical  Beach 

The  time  simulations  of  the  radiation  problem  are  conducted  until  the  force  time  history 
of  the  body  reaches  limit  state.  The  method  chosen  to  minimize  the  effect  of  transient 
behavior  associated  with  the  startup  of  the  body  from  rest  increases  the  time  required  for 
the  problem  to  reach  its  limit  state.  Either  a  large  computational  domain  or  a  method  to 
eliminate  wave  reflection  from  the  boundaries  of  the  computational  domain  is  required  to 
ensure  that  the  solution  reaches  its  limit  state. 
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The  method  developed  by  Dommermuth  and  Yue  [8]  and  used  by  Liu  [22]  for  enabling 
long  time  simulations  for  transient  ship  motions  was  used  in  this  thesis.  A  tapering  function 
was  applied  to  the  free  surface  elevation  and  free  surface  potential  to  smoothly  truncate 
these  quantities  at  the  ends  of  the  computational  domain.  The  tapering  function  is  shown 
in  Figure  5-2. 
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L 

Figure  5-2:  Representative  tapering  function. 


The  tapering  function  is  defined  as  [22]: 


1,  A  <  x  <  Z,  -  A 

Q(x,A)  =  {   n((A-ar)/A),  0  <  x  <  A 

Il((x-L  +  A)/A),  L-A<x<L 
Where  U(s)  is  the  Hermitian  polynomial  [22]  : 


(5.5) 
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U{s)  =  1  -  462s6  +  1980s7  -  3465s8  +  3080s9  -  1386s10  +  252s11    ,    0  <  s  <  1       (5.6) 

The  width  of  the  tapering  zone,  A,  is  determined  by  physical  and  computational  con- 
siderations. The  tapering  zone  must  be  wider  than  the  wavelength  of  the  longest  wave 
generated  by  the  oscillating  cylinder.  If  the  tapering  zone  is  shorter  than  this  wavelength, 
the  wave  will  be  too  strongly  "reflected"  from  the  domain  boundaries.  The  difference  be- 
tween the  length  of  the  computational  domain  and  twice  the  length  of  the  tapering  zone, 
L*  =  L  —  2 A,  must  be  long  enough  to  allow  the  simulation  to  reach  and  sustain  physically 
relevant  limit  state  behavior.  Specifically,  the  ratio  of  the  wave  length  of  the  wave  gener- 
ated by  the  moving  cylinder,  A  to  L*  must  be  small  enough  to  allow  multiple  waves  to  exist 
within  the  domain.  The  wave  length  of  the  wave  generated  by  the  forced  oscillatory  motion 
at  frequency  u  is  estimated  by  the  linear  dispersion  relation  for  gravity  waves  in  a  fluid  of 
infinite  depth,  A  =  2ng/u}2. 

Results 

The  solutions  for  the  free  surface  profile,  the  force  time  history,  and  the  amplitudes  of 
the  steady  and  oscillatory  forces  generated  by  the  combined  HOSM-BIEM  for  the  heaving 
circular  cylinder  are  compared  with  the  results  generated  from  the  other  solution  methods 
discussed  in  Chapter  1. 
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Figure  5-3  and  Figure  5-4  show  the  free  surface  profile  and  force  time  history  generated 
by  Silva  and  Peregrine  for  a  heaving  circular  cylinder.  The  parameters  for  the  problem  were 
(refer  to  Figure  5-1)  : 

•  R  =  0.25; 

•  H  =  -0.5; 

•  d  =  1.0; 


•  £  =  0.1;  and, 


•  Frequency  of  oscillation,  u;  =  3.13  rad/sec. 
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Figure  5-3:  Free  surface  profile  above  heaving  cylinder  after  approximately  5  periods  of 
oscillation.  The  ordinate  is  the  free  surface  elevation  magnified  by  a  factor  of  500.  The 
mean  position  of  the  cylinder  center  is  (0.0,  —0.5).  (Silva  and  Peregrine,  "Engineering 
Analysis  with  Boundary  Elements",  1990,  Vol.  7,  No.  4). 
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Figure  5-4:  Time  history  of  the  vertical  component  of  the  force  for  the  heaving  cylinder. 
The  ordinate  is  the  normalized  vertical  force.  The  abscissa  is  time.  (Silva  and  Peregrine, 
"Engineering  Analysis  with  Boundary  Elements",  1990,  Vol.  7,  No.  4). 
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Figure  5-5  and  Figure  5-6  show  the  free  surface  profile  and  force  time  history  gener- 
ated by  the  combined  HOSM-BIEM  method  for  a  heaving  circular  cylinder.  The  problem 
parameters  were  those  listed  on  Page  83. 
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Figure  5-5:  Free  surface  profile  above  heaving  cylinder  after  approximately  5  periods  of 
oscillation.  The  ordinate  is  the  free  surface  elevation  magnified  by  a  factor  of  500.  The 
mean  position  of  the  cylinder  center  is  (L/2,  —0.5).  (Combined  HOSM-BIEM  method). 
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Figure  5-6:  Time  history  of  the  vertical  component  of  the  force  for  the  heaving  cylinder. 
The  ordinate  is  the  normalized  vertical  force  amplitude.  The  abscissa  is  time.  (Combined 
HOSM-BIEM  method). 
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The  vertical  scale  for  the  free  surface  profiles  presented  in  Figure  5-3  and  Figure  5-5 
are  magnified  by  a  factor  of  500  to  facilitate  the  comparison  of  the  results  generated  by  the 
two  different  methods.  The  results  generated  by  the  BIEM  method  of  Pergrine  and  Silva 
and  the  results  generated  by  the  present  method  are  qualitatively  similar.  A  harmonic 
analysis  of  the  results  of  Peregrine  and  Silva  was  not  performed.  Therefore,  a  more  detailed 
comparison  of  the  results  from  the  two  methods  could  not  be  made. 

A  comparison  of  the  steady  and  first  and  second  harmonics  of  the  forces  generated  by 
the  solution  method  of  Wu  [52]  and  the  present  method  are  presented  in  Figures  5-7,  5-8, 
and  5-9.  Wu  satisfied  the  body  boundary  condition  on  the  exact  position  of  the  body  and 
used  linearized  free  surface  dynamic  and  kinematic  boundary  conditions.  Wu  expressed  the 
velocity  potential  in  terms  of  a  multipole  expansion  [52]: 


<h  =  E  E  A>r 


m=\  s=-oo 


AmO+isuit 


+  /mM,0 


+  E  E  B>r 


m—l  s—-oo 


Jm6+isujt 


+  9m(r,0>t) 


(5.7) 


eim0+iswt  jTm  jg  f^e  expression  for  the  multipole  expansion  for  a  cylinder  in  an  infinite  fluid. 
/£,  (r,  9,  t)  and  g^  (r,  9,  t)  are  required  to  enable  the  potential  to  satisfy  the  free  surface 
boundary  conditions  [52]. 
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Figure  5-7:  Vertical  drift  force  on  the  heaving  cylinder  for  fcil  =  0.1  ( )  and  kR  =  1.0 

(— )•  Wu  (1993)  (□)  ;  present  method  (Q).  (H  =  3R). 
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Figure  5-8:  First  harmonic  of  the  vertical  force  on  the  heaving  cylinder  for  kR  =  0.1  ( ) 
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Figure  5-9:   Second  harmonic  of  the  vertical  force  on  the  heaving  cylinder  for  kR  =  0.1 
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The  results  generated  by  the  method  of  Wu  and  the  present  method  agree  well  for 
heave  amplitude  to  body  radius  ratios  (£//?)  less  than  1.0.  In  the  limit  of  vanishing  heave 
amplitude  to  body  radius  ratios  both  methods  converge  to  the  linear  solution  of  Ogilve  [35]. 

The  assumptions  used  by  Wu  to  justify  using  the  linear  free  surface  boundary  condition 
become  invalid  as  the  body  approaches  the  free  surface.  For  a  deeply  submerged  body,  the 
influence  of  the  body  on  the  free  surface  is  of  order  0(R/H)  [52].  The  linear  free  surface 
boundary  condition  ignores  terms  that  involve  products  of  order  0(R/H)2  [52].  For  heave 
amplitudes  greater  than  1.0  the  body  is  close  enough  to  the  free  surface  to  invalidate  the 
deep  submergence  requirement  for  Wu's  method  to  be  valid.  In  fact,  for  the  real  (physical) 
situation  the  waves  break  for  the  large  orbital  radii  cases. 

5.2     Radiation  Problem  for  a  Circular  Cylinder  in  a  Circular 
Orbit 

The  forces  and  free  surface  profiles  resulting  from  the  interaction  of  a  fully  submerged 
cylinder  forced  to  undergo  circular  orbital  motion  with  an  initially  quiescent  free  surface 
are  presented  in  this  section.  The  forces  are  evaluated  as  functions  of: 

•  The  nondimensional  orbital  radius,  j/R; 

•  The  submergence  ratio,  h  =  H/R;  and, 

•  The  body  radius  to  wavelength  ratio  (x  2ir),  kR  —  2irR/\. 
The  forces  analyzed  are: 

•  The  first  and  second  harmonics  of  the  vertical  force,  and 

•  The  steady  vertical  force. 
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The  method  for  computing  the  forces  on  the  body  and  the  method  for  starting  the 
problem  from  rest  are  those  presented  in  the  previous  sections. 

Results 

The  free  surface  profile,  the  force  time  history,  and  the  amplitudes  of  the  steady  and 
oscillatory  forces  associated  with  the  orbiting  circular  cylinder  are  compared  with  the  results 
generated  from  the  other  solution  methods  discussed  in  Chapter  1. 

Figure  5-10  and  Figure  5-11  show  the  free  surface  profile  and  force  time  history  generated 
by  Peregrine  and  Silva  for  the  orbiting  circular  cylinder.  The  parameters  for  the  problem 
were  (refer  to  Figure  5-1)  : 

•  i?  =  0.3; 

•  H  =  -0.5; 

•  d=  1.0; 

•  Orbital  radius  =  0.025;  and, 

•  Frequency  of  oscillation,  uj  =  3.13  rad/sec. 


91 


i*<1 


0.M- 


100 


-4.M- 


<.W- 


Figure  5-10:  Free  surface  profile  above  the  orbiting  cylinder  after  approximately  5  periods 
of  rotation.  The  ordinate  is  the  free  surface  elevation  magnified  by  a  factor  of  500.  The 
cylinder  center  is  orbiting  about  the  point  (0.0,  —0.5).  (Silva  and  Peregrine,  "Engineering 
Analysis  with  Boundary  Elements",  1990,  Vol.  7,  No.  4). 
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Figure  5-11:  Time  history  of  the  vertical  ( )  and  horizontal  ( — )  components  of  the 

force  for  the  orbiting  cylinder.  The  ordinate  is  the  normalized  force  amplitude.  The  abscissa 
is  time.  (Silva  and  Peregrine,  "Engineering  Analysis  with  Boundary  Elements",  1990,  Vol. 
7,  No.  4). 
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Figure  5-12  and  Figure  5-13  show  the  free  surface  profile  and  force  time  history  gen- 
erated by  the  new  combined  HOSM-BIEM  method  for  the  orbiting  circular  cylinder.  The 
parameters  for  the  problem  were  those  listed  on  Page  91: 
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Figure  5-12:  Free  surface  profile  above  the  orbiting  cylinder  after  approximately  4  periods 
of  oscillation.  The  ordinate  is  the  free  surface  elevation  magnified  by  a  factor  of  500.  The 
cylinder  center  is  orbiting  about  the  point  (L/2  ,  —0.5).  (Combined  HOSM-BIEM  method). 
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Figure  5-13:  Time  history  of  the  vertical  ( )  and  horizontal  ( — )  components  of  the 

force  for  the  orbiting  cylinder.  The  ordinate  is  the  normalized  force  amplitude.  The  abscissa 
is  time.  (Combined  HOSM-BIEM  method). 
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As  was  the  case  for  the  heaving  cylinder,  the  comparison  with  the  results  of  Peregrine 
and  Silva  provides  a  qualitative  comparison  of  the  results  generated  by  these  two  different 
methods. 

It  has  been  shown  that  at  first  order  (c.f.  Ogilvie  [35])  and  at  second  order  (Wu  [53]) 
an  orbiting  circular  cylinder,  in  a  fluid  of  infinite  extent,  generates  waves  in  one  direction; 
e.g.,  a  cylinder  undergoing  a  clockwise  orbit  generates  only  right  going  waves.  Wu  [52] 
demonstrated  that  if  terms  up  to  fifth  order  are  included  in  a  multipole  expansion  rep- 
resentation of  the  potential,  the  possiblity  of  waves  being  transmitted  in  both  directions 
exist.  However,  these  waves  are  small.  The  results  of  Silva  and  Peregrine  and  the  present 
results  demonstrate  that  in  shallow  water  waves  of  significant  amplitude  are  generated  and 
transmitted  in  the  "upstream"  direction. 

This  result  is  in  contrast  to  that  for  a  cylinder  undergoing  forced  orbital  motion  in  a 
fluid  of  infinite  depth.  Figure  5-14  shows  the  free  surface  profile  above  a  circular  cylinder 
orbiting  beneath  the  free  surface  in  a  fluid  of  infinite  depth.  The  problem  parameters,  with 
the  exception  of  the  fluid  depth  were  those  listed  on  Page  91.  The  cylinder  generates  waves 
in  only  the  "downstream"  direction  when  the  fluid  depth  is  infinite. 
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Figure  5-14:  Free  surface  profile  above  orbiting  cylinder  in  a  fluid  of  infinite  depth  after 
approximately  5  periods  of  oscillation.  The  ordinate  is  the  free  surface  elevation  magnified 
by  a  factor  of  500.  The  cylinder  is  orbiting  about  the  point  (L/2  ,  —0.5).  (Combined 
HOSM-BIEM  method). 
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A  comparison  of  the  steady  and  first  and  second  harmonics  of  the  forces  generated  by 
the  solution  method  of  Wu  [52]  and  the  present  method  are  presented  in  Figures  5-15,  5-16, 
and  5-17. 
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Figure  5-15:  Vertical  ( )  and  horizontal  ( — )  drift  forces  on  orbiting  circular  cylinder. 

Wu  (1993)  (□);  present  method  (Q).  (H  =  3R,  kR  =  0.5). 


As  was  the  case  for  the  heaving  cylinder  the  two  methods  agree  well  for  small  amplitude 
of  oscillation  (orbital  radius)  to  cylinder  radius  ratios. 

For  large  orbital  radii,  the  cylinder  generates  waves  with  local  steepnesses  that  exceed 
the  limiting  Stokes  steepness.  Figure  5-18  shows  the  local  wave  steepness  for  a  case  where 
the  orbital  radius  to  body  radius  is  1.5  and  the  mean  depth  is  3  times  the  body  radius.  For 
large  orbital  radii,  the  method  becomes  unstable  as  the  body  motion  generates  excessively 
steep  waves,  and  does  not  reach  limit  state  behavior  before  the  simulation  stops. 
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Figure  5-16:  First  harmonic  of  the  vertical  ( )  and  horizontal  ( — )  forces  on  orbiting 

circular  cylinder.  Wu  (1993)  (D);  present  method  (Q).  (H  =  3R,  kR  =  0.5). 


Figure  5-17:  Second  harmonic  of  the  vertical  ( )  and  horizontal  ( — )  forces  on  orbiting 

circular  cylinder.  Wu  (1993)  (□);  present  method  (Q).  (H  =  3R,  kR  =  0.5). 
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Figure  5-18:  Free  surface  profile  above  the  orbiting  cylinder.  (H  =  3R;  kR  =  0.5). 
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5.3     Summary 

The  ability  of  the  combined  high-order  spectral-  boundary  integral  equation  method  to  sim- 
ulate large  amplitude  oscillatory  motions  of  a  circular  cylinder  beneath  a  free  surface  has 
been  demonstrated.  The  method  enables  forces  and  free  surface  profiles  to  be  evaluated  up 
to  limiting  Stokes'  wave  steepness  for  a  range  of  problems  of  interest.  Problems  involving 
either  finite  or  infinite  depth  can  be  modelled.  The  importance  of  using  a  solution  method 
which  identifies  the  existence  of  physically  impossible  waves  is  demonstrated  by  the  applica- 
tion of  the  present  method  to  large  amplitude  body  motions  near  a  free  surface.  A  method 
such  as  the  one  used  by  Wu,  which  does  not  account  for  non-linear  free  surface  interactions 
and  locally  steep  waves,  extends  linear  results  into  flow  regimes  where  linear  assumptions 
are  invalid. 

In  the  next  chapter,  the  HOSM-BIEM  method  will  be  applied  to  the  combined  radiation- 
diffraction  problem  for  a  fully  submerged  cylinder  free  to  respond  to  incident  waves. 
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Chapter  6 


Combined  Radiation  and 


Diffraction  Problem  for  a 


Submerged  Circular  Cylinder 


The  purpose  of  this  chapter  is  to  demonstrate  the  application  of  the  combined  high  order 
spectral-boundary  integral  equation  method  to  the  problem  of  a  fully  submerged  circular 
cylinder  free  to  respond  to  an  incident  wave  field.  Figure  6-1  defines  the  problem  parameters. 
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Figure  6-1:  Problem  parameters  for  the  combined  problem 


•  R  =  the  body's  radius 


r0  =  the  radius  of  body's  orbital  trajectory 


(x0,  h0)  =  the  location  of  the  point  about  which  the  body  is  orbitting 


A  =  the  wavelength  of  the  ambient  waves 


A  =  the  wave  amplitude 
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6.1     The  Combined  Radiation  and  Diffraction  Problem  for  a 
Circular  Cylinder  Free  to  Respond  to  Waves 

The  analyses  of  neutrally  bouyant  circular  cylinders  free  to  respond  to  incident  wave  fields 
were  facilitated  by  the  generation  of  cylinder  trajectories.  The  cylinder  trajectories  were 
generated  by: 

•  Specifying  initial  conditions  for  the  cylinder  position,  velocity  and  acceleration; 

•  Solving  for  the  instantaneous  force  on  the  cylinder  using  the  solution  methodology 
detailed  in  Chapters  2  and  3;  and, 

•  Updating  the  cylinder  position,  velocity,  and  acceleration  using  the  state  varible  ap- 
proach represented  by  Equations  2.33,  2.34,  2.35,  2.36,  and  2.37. 

Initial  Conditions 

Ogilvie  [35]  used  linear  theory  to  show  that  the  amplitude  and  phase  angle  of  the  orbital 
motion  of  a  submerged  circular  cylinder  differs  from  that  of  a  water  particle  centered  at 
the  location  of  the  cylinder's  origin  by  a  quantity  of  fourth  order  in  kR,  where  k  is  wave 
number  and  R  is  the  body's  radius.  This  fact  was  used  to  estimate  an  intitial  velocity  and 
accleration  for  the  cylinder. 

The  first-order  orbital  radius  of  a  water  particle  located  initially  at  (x0,  h0)  (see  Figure  6- 
1)  would  be: 

ro(0)  =  Ae-kh°e-twt\t=o  =  Ae~kh°  (6.1) 

where  u  is  the  wave  frequency. 

104 


The  initial  velocity  (v0{t  =  0)  =  (u(0),iu(0)))  and  acceleration  {a0(t  =  0)  =  (oz(0),ay(0))) 
of  the  cylinder  were  set  the  same  as  that  of  a  fluid  particle  at  the  cylinder  center. 


u(0)  =  uAe-kh°  cos(kx)  (6.2) 

u;(0)  =  u)Ae~kh°  sin(fca:)  (6.3) 

ax(0)  =  uj2Ae-kh°  sin(fcx)  (6.4) 

ay(0)  =  -J2 Ae~kho  cos(hx)  (6.5) 

Method  of  Computing  Forces  and  Generating  the  Trajectory  of  the  Cylinder 

The  potential  on  the  body  is  computed  at  each  time  step  by  solving  the  initial-boundary 
value  problem  presented  in  Chapters  2  and  3.  The  potential  is  used  to  compute  the  instan- 
taneous force  on  the  body,  F(t),  by  the  application  of  Equation  5.3.  The  computed  force 
is  used  to  update  the  body's  acceleration,  velocity,  and  position  using  the  state  variable 
representation  of  these  quantities  presented  in  Chapter  2: 


a0{t)    =    X(t)  =  F{t)/M  (6.6) 

rt+6t 
v0{t  +  5t)     =    v0(t)  +  /        a0{t)  dt  (6.7) 

rt+St 

a0(t  +  St)     =    a0(t)+  v0(t)dt  (6.8) 

Numerical  Results 

The  trajectories  of  submerged  cylinders  beneath  waves  were  simulated.  The  waves  used 
in  the  simulations  were  waves  travelling  in  the  positive  X  direction  (left  to  right  in  the 
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figures  that  follow). 

The  convergence  with  respect  to  time  step  size  on  the  trajectory  of  a  cylinder  free  to 
respond  to  waves  is  shown  in  Figure  6-2.  Figure  6-2  shows  that  in  order  to  achieve  accurate 
cylinder  trajectories  a  time  step  size  of  T/256  or  smaller  is  required. 
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Figure  6-2:  Convergence  of  the  trajectory  of  a  submerged  circular  cylinder  with  respect  to 
time  step  size,  dt.  Wave  slope  kA  =  0.08;  cylinder  initial  position  (x0,h0)  =  (157T,  —1.2); 
kR  =  0.4.  The  trajectory  of  the  cylinder  center  is  shown.  Order  M  =  4,  number  of  Fourier 
modes  Nx  =  2048,  domain  size  Nw  =  16,  Nfuter  =  24Nw,  Nb  =  128. 
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The  convergence  with  respect  to  number  of  body  segments  on  the  trajectory  of  a  cylinder 
free  to  respond  to  waves  is  shown  in  Figure  6-3. 
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Figure  6-3:  Convergence  of  the  trajectory  of  a  submerged  circular  cylinder  with  respect  to 
number  of  body  segments,  -/V#.  Wave  slope  kA  =  0.08;  cylinder  initial  position  (x0,h0)  = 
(157r,  — 1.2);  kR  =  0.4.  The  trajectory  of  the  cylinder  center  is  shown.  Order  M  =  4, 
number  of  Fourier  modes  A^  =  2048,  domain  size  Nw  =  16,  Nfuter  =  24iVV,  dt  =  T/256. 
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Figures  6-4  and  6-6  show  the  trajectories  for  two  different  cylinders  free  to  respond  to 
incident  waves.  The  body  radius  for  the  problem  represented  by  Figure  6-4  is  2.5  times  the 
body  radius  for  the  problem  represented  by  Figure  6-6.  The  expanded  cylnder  trajectories 
for  each  case  are  shown  in  Figures  6-5  and  6-7.  Figures  6-4  and  6-6  show  the  cylinder 
positions  at  three  different  time  steps.  The  free  surface  profiles  shown  are  the  free  surface 
profiles  above  the  cylinder  when  the  cylinders'  positions  are  as  shown  by  the  solid  lines. 
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Figure  6-4:  Cylinder  trajectory  and  an  associated  free  surface  profile.  Cylinder  position  at 

t  =  0  ( );  cylinder  position  at  t  =  3T  ( — );  cylinder  position  at  t  =  3.75T  (—•  —  ■—). 

Wave  slope  kA  =  0.08;  kR  =  1.0;  cylinder  initial  position  (x0,h0)  =  (L/2,  — 1.25).  Order 
M  —  4,  number  of  Fourier  modes  Nx  =  2048,  domain  size  Nw  =  16,  Nfuter  =  24./VV, 
dt  =  T/256,  NB  =  192. 
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Figure  6-5:  Expanded  cylinder  trajectory.  Wave  slope  kA  =  0.08;  kR  =  1.0;  cylinder  initial 
position  {x0,h0)  =  (L/2,  — 1.25).  Order  M  =  4,  number  of  Fourier  modes  Nx  =  2048, 
domain  size  Nw  =  16,  Nfilter  =  247VV,  dt  =  T/256,  NB  =  192. 
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Figure  6-6:  Cylinder  trajectory  and  an  associated  free  surface  profile.  Cylinder  position  at 

t  =  0  ( );  cylinder  position  at  t  =  3T  (—  •  —  ■  — );  cylinder  position  at  t  =  IT.  Wave 

slope  kA  =  0.08;  kR  =  0.4;  cylinder  initial  position  (x0,  h0)  =  (L/2,  —1.25).  Order  M  =  4, 
number  of  Fourier  modes  -/Vx  =  2048,  domain  size  Nw  =  16,  Njuter  =  2ANw,  dt  =  T/256, 
NB  =  192. 
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Figure  6-7:  Expanded  cylinder  trajectory.  Wave  slope  kA  =  0.08;  kR  =  0.4;  cylinder  initial 
position  (x0,  h0)  =  (L/2,  —0.8).  Order  M  =  4,  number  of  Fourier  modes  Nx  =  2048,  domain 
size  Nw  =  16,  Nfilter  =  24NW,  dt  =  T/256,  NB  =  128. 
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The  cylinder  trajectories  reflect  the  influence  of  the  imposed  initial  conditions  in  that 
for  both  cases  the  cylinders  initially  follow  the  trajectory  that  would  be  assumed  by  a  fluid 
particle  located  at  the  cylinder's  initial  position  (x0,y0).  As  time  increases  the  steady  and 
oscillatory  forces  induced  by  the  fluid-body  interaction  generate  motions  which  differ  from 
that  of  a  fluid  particle.  For  the  large  body  the  presence  of  a  negative  horizontal  drift  motion 
is  apparent  while  for  the  smaller  body  this  does  not  appear  to  be  the  case. 

The  simulation  of  the  interaction  between  the  large  cylinder  and  the  free  surface  was 
terminated  due  to  the  development  of  short,  locally  steep  waves  as  the  body  became  close 
to  the  free  surface.  The  presence  of  the  small  locally  steep  waves  is  apparent  in  Figure  6-4. 
The  development  of  locally  steep  waves  helps  to  define  the  limitations  of  the  method  of 
this  thesis  and  to  identify  flow  regimes  where  linear  theories  used  for  modelling  wave-body 
interactions  are  invalid.  The  free  surface  boundary  conditions  used  in  this  thesis  can  not 
be  used  to  simulate  breaking  waves. 

6.2     Summary 

The  combined  high-order  spectral-boundary  integral  equation  method  has  been  successfully 
applied  to  the  problem  of  a  fully  submerged  neutrally  bouyant  cylinder  free  to  respond  to 
incident  waves.  The  method  is  stable  enough  to  enable  simulations  of  sufficient  duration  to 
evaluate  the  steady  and  oscillatory  responses  of  interest. 
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Chapter  7 


Conclusions  and  Recommendations 


A  combined  high-order  spectral  method  boundary  integral  equation  method  was  developed 
for  the  analysis  of  the  nonlinear  interaction  of  non-breaking  waves  with  fully  submerged 
bodies  of  arbitrary  geometry  undergoing  arbitrary  motions.  The  motivation  for  the  devel- 
opment of  the  method  was  the  desire  to  improve  the  computation  speed  for  wave  interaction 
problems  involving  bodies  of  arbitrary  geometry  undergoing  arbitrary  motions  near  the  free 
surface. 

The  method  developed  in  this  thesis  represents  the  total  fluid  potential  as  an  expansion 
in  a  perturbation  series  of  terms  proportional  to  wave  slope.  The  total  potential  is  the  sum 
of  a  spectral  potential  and  a  body  potential.  The  method  uses  a  spectral  representation  of 
free  surface  quantities  and  a  boundary  element  representation  of  the  body.  The  spectral 
representation  of  the  free  surface  elevation  and  potential  enables  the  fast  computation  of 
free  surface  quantities.  The  boundary  element  representation  of  the  body  enables  bodies  of 
arbitrary  shape  to  be  modelled.  Fully  nonlinear  free  surface  boundary  conditions  are  used 
as  time  evolution  equations  for  the  free  surface  elevation  and  free  surface  potential. 

The  effectiveness  of  the  method  was  demonstrated  through  the  study  of  three  classes 
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of  problems:  diffraction  of  incident  waves  by  a  stationary  body,  radiation  from  a  body 
undergoing  forced  oscillatory  motion  beneath  a  free  surface,  and  a  combined  radiation- 
diffraction  problem  for  a  body  free  to  respond  to  an  incident  wave  field. 

The  success  and  limitations  of  the  present  method  in  computing  the  forces  associated 
with  nonlinear  wave  body  interactions  are  demonstrated  in  Chapters  4,  5  ,  and  6.  The 
present  method  succeeds  in  computing 

•  Oscillatory  forces; 

•  Steady  drift  forces; 

•  Free  surface  profiles  up  to  limiting  wave  steepness;  and, 

•  Trajectories  for  bodies  free  to  respond  to  waves. 

The  forces  and  free  surface  profiles  produced  by  the  present  method  for  the  diffraction 
and  radiation  problems  compare  well  with  those  produced  by  other  methods  restricted  to 
interactions  involving  non-breaking  waves,  and,  for  waves  of  small  steepness,  compare  well 
with  those  produced  by  linear  methods.  Through  evaluation  of  the  free  surface  profiles 
generated  by  the  wave-body  interaction,  locally  steep  profiles  which  exceed  limiting  wave 
steepness  can  be  identified  and  used  to  determine  when  methods  based  on  perturbation 
expansions  are  no  longer  appropriate  or  accurate  for  the  characterization  the  posed  physical 
problems. 

The  ability  of  the  present  method  to  evaluate  bodies  of  arbitrary  geometry  was  demon- 
strated by  the  evaluation  of  a  body  of  pontoon  cross  section.  The  ability  of  the  present 
method  to  evaluate  bodies  undergoing  arbitrary  oscillatory  motions  was  demonstrated  by 
evaluating  the  large  amplitude  oscillatory  motions  of  a  circular  cylinder. 
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Recommendations 

The  present  method  is  formulated  to  accomodate  non-breaking  waves.  As  shown  in  the 
analysis  of  the  diffraction  and  radiation  problems,  locally  steep  waves  limit  the  ability  of 
the  present  method  to  evaluate  wave-body  interactions  when  the  body  is  close  to  the  free 
surface  even  when  the  waves  in  the  ambient  wave  field  are  not  steep. 

The  periodic  nature  of  the  formulation  of  the  present  method  limits  problem  simulation 
time.  For  the  radiation  problem  the  use  of  a  numerical  damping  zone  enables  longer  simlu- 
ation  times.  A  method  for  enabling  long  time  simulations  for  the  diffraction  and  combined 
problems  needs  to  be  developed  to  enable  the  method  to  perform  long  time  simulations. 
This  is  important  for  the  manuevering  simulations  where  the  distance  travelled  by  a  body 
is  more  than  a  few  body  lengths. 

The  present  method  was  demonstrated  for  two  dimensional  non-lifting  bodies.  The 
extension  of  the  method  to  three  dimensions  and  problems  involving  lift  is  in  principle 
straightforward,  but  remains  to  accomplished  and  demonstrated. 
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